ON THE ONE-DIMENSIONAL STABILITY OF VISCOUS 
STRONG DETONATION WAVES 
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Abstract. Building on Evans function techniques developed to study 
the stability of viscous shocks, we examine the stability of viscous strong 
detonation wave solutions of the reacting Navier-Stokes equations. The 
primary result, following ^ I17| . is the calculation of a stability index 
whose sign determines a necessary condition for spectral stability. We 
show that for an ideal gas this index can be evaluated in the ZND limit 
of vanishing dissipative effects. Moreover, when the heat of reaction is 
sufficiently small, we prove that strong detonations are spectrally stable 
provided the underlying shock is stable. Finally, for completeness, the 
stability index calculations for the nonreacting Navier-Stokes equations 
are included. 
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1. Introduction and Preliminaries 

1.1. Introduction. Laboratory and numerical experiments indicate that 
detonations have quite sensitive stability properties. Indeed steady planar 
detonations subjected to one-dimensional longitudinal perturbations may 
change form to "galloping" detonations in which the velocity fluctuates pe- 
riodically in time. Such detonations have been been predicted numerically 
|15j and observed experimentally in various settings by [22], and |50j . 
Another instability, this one with 3-dimensional structure, is the "spinning 
detonation" long- known in lab experiments El an d more recently captured 
numerically for the ZND model in both [3] and |28| . A three-dimensional 
perturbation of a steady detonation wave propagating down a tube with 
a circular cross-section may bifurcate to a wave with a complex rotating 
structure which traces a helical path along the boundary of the tube. This 
structure is typically followed by localized regions of extremely high pres- 
sure. 

Due to this sensitivity and the complicated, nonlinear nature of react- 
ing gas dynamics (which includes such difficulties as turbulence, boundary 
layers, and complex chemical interactions), stability analyses of detonation 
waves are largely numerical studies of the ZND model (which neglects dis- 
sipative effects) as in e.g. , [SB > EH1 E3 5 or are restricted to various 
incarnations of the Majda or Majda-Rosales 1 models (unphysical analogues 
of Burgers equation) as in (S2J , jM] , [HE] , [HSI EH1 EEl EH EZ] , ,and [EE]. Our 
approach, utilizing the Evans function, allows the treatment of the reacting 
Navier-Stokes equations and yields an explicitly computable quantity known 
as the stability index. 

In Eulerian coordinates the reacting Navier-Stokes equations modeling 
the simplest possible one-step chemical reaction can be written 

(1.1) Pt + (pu) x = 0, 

(1.2) (pu) t + (pu 2 +p) x = (uu x ) x , 

(1.3) £f + [puS + up] x = {6T X ) X + (qpdY x ) x + {vuu x ) x , 

(1.4) (pY) t + {puY) x = (pdY x ) x -k P Yv(T). 

In (|l.lj) - (|1.4|) and below we use, unless stated otherwise, the notations: p, 
u, p, £, T, and Y represent respectively density, velocity, pressure, total 
energy, temperature, and mass fraction of reactant. The use of the tilde 
denotes that the energy £ = p(u 2 /2 + e) is modified from the standard gas- 
dynamical energy £ = p(u 2 /2+e) = pE due to heat produced in the chemical 
reaction by e = e + qY . The positive constants v, 9, and d represent the 
effects of viscosity, heat conductivity, and species diffusion. The positive 
constants k and q measure the rate of reaction and the heat released in 
reaction, respectively, and the form of the so-called ignition function (p is 



^For the remainder of the paper, we refer to all models with the simplifying feature of 
scalar kinetics as the Majda model. 
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discussed in detail below. The system is closed by specifying equations of 
state, p = p(p, e, Y) and T = T{p, e, Y). We begin by assuming only that p 
and T are independent of Y, but for some portions of the analysis we shall 
assume further an ideal, polytropic gas, i.e 

P{P, e) = Tpe, T(p, e) = c~ l e, 

where c v , the specific heat at constant volume, and T, known as the Gruneisen 
coefficient, are constants. Equations (|l.ljl - (|1.4jl are standard; a derivation 
can be found in |UT] . 

Often ip is assumed to satisfy the Arrhenius law, so 

V?(T) = exp 

where Ea is the activation energy and R is the gas constant (assuming the 
ideal gas law, R = c v T ). However, nonvanishing of the exponential creates a 
problem known as the "cold-boundary difficulty." Essentially, nonvanishing 
of ip precludes the unburned state from being a rest point of the traveling 
wave ODE. In place of the Arrhenius kinetics we make the standard as- 
sumption that the smooth function <p vanishes for temperatures below some 
ignition temperature, Tj, and is identically 1 for some larger value of T. This 
circumvents the cold boundary difficulty; see Figure lTTl This model includes 




T i T 



Figure 1.1. The Ignition Function ip 

dissipative effects neglected by the ZND model, and allows for complete gas 
dynamical effects unlike the Majda model. We remark that an artificially 
strictly parabolic multidimensional version is considered in the appendix of 
|65| : here we include the additional difficulty of partial parabolicity. 

1.2. Plan of the Paper. In Section^we provide the relevant background 
material for our analysis. We first discuss the standard ZND model as a 
prelude to a discussion of the analysis of JH]- ZND solutions are singu- 
lar solutions in the context of geometric singular perturbation theory used 
therein. The structure of these singular solutions allows us to evaluate the 
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stability index, which we do in Sectional The backgound material in Sec- 
tion H concludes with a description of the Evans function theory for the 
stability of viscous shock waves. Section |2] contains the stability index cal- 
culations for the (nonreacting) Navier-Stokes equations. These prove useful 
in Section 12 where we examine the reacting system (|1.1[) - ([1.4|) . In Appen- 
dix El we include a revised version of an appendix of which is used in 
the stability index calculations. 

1.3. The ZND Model. Setting the constants v, 9, and d equal to zero in 
(jl.ljl - p.4jl yields the ZND model introduced independently by Zeldovich, 
von Neumann, and Doring. In their formulation dissipative effects are ne- 
glected and the reaction rate is assumed to be finite. This is a refinement 
of the early Chapman-Jouget model in which the reaction was assumed to 
take place instantaneously. The ZND model, in Eulerian coordinates, then 
has the form 

Pt + (pu) x = 0, 
(pu) t + (pu 2 + p) x = 0, 
£ t + ({£ + P )u) x = 0, 
( P Y) t + {puY) x = -kpY<p(T). 

As above the constant k is the reaction rate, and </? is the ignition function. 
Strong detonations in the ZND model are initiated by a (purely) gas dynam- 
ical shock, called the Neumann shock, which heats the gas by compressing 
it. The increase in the temperature to a sufficiently high level "turns on" ip 
and thus starts the reaction. Thus these waves have the structure of a gas 
dynamical shock followed by a reaction zone resolving to the final burned 
state. This is seen in the characteristic "detonation spikes" in the temper- 
ature and pressure profiles in strong agreement with observed features in 
laboratory experiments. We'll see this structure in our discussion of j!8j 
below. 

1.4. Existence of Strong Detonations. In j!8j . existence of traveling 
wave solutions of Equations (jl.l|) - (jl.4j) was studied using the techniques 
of geometric singular perturbation theory (GSPT). As the orientations of 
the singular manifolds constructed in that argument will play a role in our 
analysis, we recap the argument here. An interesting feature of the GSPT 
analysis of (jl.l[) - (jl.4j) is the recovery of the shock layer analysis of [2T] . 

1.4.1. The Hugoniot Curve. Traveling wave solutions of Equations p.lj) - 
p.4j) are those which depend only on £ = x — st. This ansatz reduces the 
system p.ljl - p.4jl to a system of ordinary differential equations (ODEs). 
By Galilean invariance we may, without loss of generality, set s = 0, so the 
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system is 

(1.5) (pu)' = 0, 

(1-6) (pu 2 +p)' = (uu'Y, 

(1.7) fpu(y + e)+up) = (0T')' + (qpdY')' + (vuu')\ 

(1.8) (puF)' = (pc/YO'-fyrtVCT), 

where ' denotes differentiation with respect to x. From the first equation 
(|1.5|) it follows that the mass flux m = pu has a constant value. Moreover 
each of ()1.6|) and (|1.7|) can be integrated up once. We suppose that an 
unburned state p + ,u + ,p + ,Y + has been fixed at +oo. Then the momentum 
equation integrates to 

pu 2 + p — {p+u\ + p + ) = vu'. 

For a possible connection to a burned state p, u,p, Y at — oo, it is necessary 
that the state be a rest point of the ODE, or more precisely p,u,p,Y must 
satisfy 

pu 2 +p- {p+u\ + p + ) = 0. 
Searching for all such states leads to the expression 

p — p + = —m I 

This equation describes a line in the specific volume-pressure plane with 
slope —m 2 . It is referred to as the Raleigh Line. Similarly integrating the 
energy equation and searching for rest points again, leads to the equation 
for the Hugoniot Curve 

~ ~ 1 , U 1 1 

e - e+ = -~{p + p+) 

2 \p P+ 

The intersection of the Raleigh line and the Hugoniot curve in the specific 
volume-pressure plane determine the possible burned states corresponding 
to the fixed unburned state p + , u + , p + , Y + . An important distinction from 
the nonreacting case is the fact that there may be one, two, or no possible 
burned states. See Figure IT~21 Another interesting feature of Figure IT~2l is 
the fact that the Hugoniot curve splits into two branches. This indicates 
that the conservation relations are compatible with two distinct types of pro- 
cesses, just as observed by early experimentalists. The compressive solutions 
are called detonations while the expansive solutions are referred to as defla- 
grations. Accordingly the two branches of the Hugoniot curve are called the 
detonation branch and deflagration branch. Each of these branches is fur- 
ther subdivided into two sections. Here we focus on the detonation branch, 
but similar characterizations apply to the deflagration branch. From the 
diagram it's clear that there is a unique value of —m 2 so that the Raleigh 
line is tangent to the Hugoniot curve. The point of tangency is called the 
Chapman-Jouget point, and a detonation connecting to the burned state 
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Figure 1.2. The Hugoniot Curve 



identified by that point is called a Chapman-Jouget detonation. For values 
of —m 2 which are smaller than this unique value, the Raleigh line intersects 
the Hugoniot curve twice. The larger of these two intersections is the burned 
state corresponding to a strong detonation while the smaller corresponds to 
the burned state of a weak detonation. One distinction between these waves 
is the following: strong detonations satisfy the Lax characteristic condition 
while weak detonations are under compressive. Here compressivity refers to 
the number of incoming characteristics. Finally, if —m 2 is too large, the 
Raleigh line and the Hugoniot curve do not intersect and there are no pos- 
sible burned states which are compatible with the unburned state. 



1.4.2. GSPT Analysis of Detonation Waves. We begin with the briefest of 
introductions to GSPT. Consider a system of singularly perturbed ODE 



(Slow e 



f = /(*,¥). 
e%=9(.x,y), 



where e is small. We call such a system, the slow system. By rescaling the 
independent variable by r = t/e, we obtain the equivalent (when e / 0) fast 
system 



(Fast e 
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Looking at flSlowJ and flFastJ , it's clear that there are then two distin- 
guished limiting systems when e = 0. They are the reduced problem 



(SIowq) 

and the layer problem 
(Fasto) 



= g(x,y), 



f = 0, 



The basic idea then is to construct solutions of the original system as smooth 
perturbations of the composite orbits of the decoupled limiting systems 
flSlowoD and flFastoD - For more details consult ^SI and |55] . 

With this framework in mind, we take up the analysis of [TS|. We note 
that here, following |18j . we are assuming an ideal, polytropic gas, so that 
p = RpT, e = c v T, and T = 7 — 1 = R/c v . We remark that existence 
for more general equations of state has been shown by different methods 
in PI], but the singular manifolds constructed in the GSPT analysis ^H] 
play a key role in the evaluation of our stability condition. Namely, they 
contain the necessary geometric information about the profile. Using m 
and the integrated versions of (jl.6|) and l|1.7|) . one finds by some further 
manipulation that 

(T X 1 -)- 

u u± 

6T x + vuu x +qpdY x = m({R + c v ){T -T±) + q(Y -Y±) + ^(u 2 - u 2 ± 

Next define the variable Z by the relationship 

Y Y 
Z = Y-pd^ = Y-d^, 
m u 

and note that Y x vanishes at ±00 to obtain 

z_ = y_, z+ = y + . 

The equation for Y can thus be rewritten as 

Y 

Z x = -k—(p(T). (Note u ^ 0) 
u 
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system 




(1.9) 


vu x = 


(1.10) 


0T X = 


(1.11) 


dY x = 


(1.12) 


Z x = 



Finally, rescaling to make the equations dimensionless, one arrives at the 



1 (T 

T - 1 [T -u)+qZ (u - 1) , 

7 2 

u(Y-Z), 
Y 
u 

Here all quantities have been rescaled, M, defined by M 2 = u 2 / '(pjRT), is 
the Mach number. 

The values of the dissipative coefficients (u,0, and d) are typically quite 
small. Taking advantage of this smallness, the next step is to fix small values 
z>, 9, and d, and then to set v = eu, 6 = e9 and d = ed so that the system 
(01>- (fTTT2|) takes the form 

(1.13) dto. = „- 1 + -^(I-l 

(1.14) e§T x = T-l- 1 —{T-u)+qZ- (7 ~ 1)M " (u - l) 2 , 

7 2 

(1.15) edV x = w(y-Z), 

(1.16) = <p(T). 

u 

Here e is supposed to be small, and this system is singularly perturbed. 
Setting e = yields the reduced (slow flow) system 

(1.18) = T-l--^(T-u)+qZ- ^~y M (u-1) 2 , 

7 2 

(1.19) = u(Y-Z), 

(1.20) Z x = —<p(T). 

u 

Equations (|1.17|) - (|1.19|) define a one-dimensional manifold C upon which 
equation (|1.20|) describes a flow. Upon noting that (jl,17j) is independent of 
Y and Z; and (|1.18|) is independent of Y; and when it / 0, (jl,19j) implies 
that Y = Z; C can be visualized in three-dimensional uTZ-space. The 
equation (jl,17j) describes a parabolic trough. Using (|1.17j) in (|1.18j) yields 

= Tl±l +q Z + u^±(l + 1 M 2 )-l-^- 1)M2 , 
27 H 27 ^ ' ' 2 
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which describes a plane K, in uTZ-sp&ce. The manifold C is exactly the 
intersection of this plane and the parabolic trough. This intersection is pic- 
tured in Figure E3 Note that C splits into two branches, as the requirement 
that there be two burned end states corresponding to the fixed unburned 
state forces the vertex of the intersection to have a negative Z coordinate. 
In fact the vertex is exactly at Z = in the boundary Chapman-Jouget 
case. Note also that all the rest points of l)1.13|) - l)1.16jl are contained in C. 




Figure 1.3. The Intersection of K, and the Trough 



Rescaling the independent variable by £ = x/e in (|1.13j) - p.l6|) yields the 
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equivalent (when e / 0) system 
(L21) ^ = + 

(1.22) 9T £ = T — 1 — 2^1(T - u )+qZ- ^~ (u - l) 2 , 

7 2 

(1.23) dY ? = u(Y-Z), 

(1.24) Z ( = e (~<p(T)Y 

Setting e = in (|1.21)) - (|1.24|) yields the layer (fast flow) system 

(1.25) 0u t = „-i + _^(I-l), 

(1.26) = T-1-- {T-u)+qZ-— — ^ (u-l) 2 , 

(1.27) = u(Y-Z), 

(1.28) Z ? = 0, 

which is the gas dynamic shock problem of [2^ coupled to (jl.27j) . Thus 
Z 

fast flow 

" slow flow 




Figure 1.4. Singular Flow 
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each constant Z slice, the fast flow is described by the shock layer analysis of 
|21j . Figure IQ1 represents the structure of the singular (e = 0) flow looking 
down on C from a vantage point perpendicular to the plane /C. Hollow arrows 
represent fast flow, while single arrows represent slow flow. The fast flow in 
the plane Z = constant corresponds to a nonreacting gas dynamical shock, 
while the slow flow proceeds on each branch of C representing the progress of 
the reaction. Combining, we see in Figure 1131 the perturbed composite orbit 
of a strong detonation. Note the presence of the ZND structure, namely the 
gas dynamical shock to the Neumann spike which raises the temperature 
above ignition followed by a reaction resolving to the final totally burned 
state. The same structure can be seen in the analysis of the scalar Majda 
model EU SH • Having dealt with the question of existence (at least 



for small u, 9, and d), we turn our attention to the tools of our stability 
analysis. 

1.5. The Evans Function and The Gap and Tracking Lemmas. The 

proper notion of stability for traveling waves connecting constant end states, 
as are the detonations and shocks we consider, is that of orbital stability, 
that is the convergence of the perturbed solution to the manifold of solu- 
tions which connect the same two end states. We refer to this as nonlinear 
stability. A weaker notion is that of linearized orbital stability defined as 
convergence of the perturbation solving the linearized equations to the tan- 
gent manifold of solutions connecting the end states. Closely related to this 



T 



unburned state 




z=0 



u 



Figure 1.5. Detonation 
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concept of linear stability is spectral stability; a profile is spectrally stable 
if the linearized operator about the wave, L, has no spectrum in the set 
{A £ C | Re A > 0} except A = 0. We note that translation invariance 
implies that is an eigenvalue of the linearized operator, so necessarily 
£ &(L). Spectral stability is clearly necessary for linear stability, itself 
necessary for nonlinear stability. Recent work by Mascia and Zumbrun 
|46| 1471 148j extending earlier work of Zumbrun and Howard [HE1 indicates 
that spectral stability implies nonlinear orbital stability in the settings of 
viscous conservation laws and relaxation systems which are closely related 
to the combustion systems we consider. In light of these results, the deter- 
mination of spectral stability can be regarded as the essential initial step in 
determining the stability of detonation waves. We also note that |621 I67| 
contain large-amplitude nonlinear stability results, relevant to the detona- 
tion problem. 

1.5.1. Background. As noted above, a vital step in determining stability 
amounts to locating the spectrum of a linear operator. The search for 
spectrum is facilitated by the Evans function, D(\), a potent tool in the 
investigation of stability of traveling waves. This function, analytic on the 
unstable half-plane, is an infinite-dimensional analogue of the characteristic 
polynomial. Zeros of D(-) correspond to eigenvalues of the linearized oper- 
ator about the wave. The Evans function was introduced in |3 1 1 0| lll| 11 2 j 
specifically to study nerve axon equations and further developed in 11 to the 
case of semilinear parabolic systems. In |17j the use of the Evans function 
was extended to the case of 2 x 2 viscous conservation laws. The extension 
tonxji systems has been completed in [2]. 

The construction of the Evans function is accomplished by analytically 
parametrizing the unstable/stable manifolds of the variable coefficient eigen- 
value equations. This is done by comparing these objects to the correspond- 
ing unstable/stable manifolds for the limiting constant coefficient systems 
at rtoo. We now give an abbreviated description of the construction for the 
case of a system of viscous conservation laws 

U t + f(U) x = (B(U)U x ) x , 

where U, f € R n and B is an n x n matrix. As we shall be concerned 
with the cases of gas dynamics and combustion, it will be the case that the 
matrix B is incompletely parabolic. However, the equations of compressible 
gas dynamics satisfy the symmetrizability, dissipativity, and block structure 
conditions of Kawashima |30| : 
symmetrizability 

There exists a symmetrizer A°(U), symmetric and positive definite, 
(+) such that A°(U)A(U) is symmetric and A°(U)B(U) 
is symmetric and positive semidefinite. 



11 



GREGORY LYNG AND KEVIN ZUMBRUN 



dissipativity 

(++) There is no eigenvector of A(U) lying in the kernel of B(U). 
block structure 

(+++) The right kernel of B(U) is independent of U. 

In the above A(U) denotes the Jacobian matrix of the flux /. We discuss 
the construction in the case of combustion in more detail below. A viscous 
profile, 

U(x,t) = U(x - st), U(±oo) = U±, 
is a solution of the (integrated) traveling wave ODE 

_ s(C7 _ + fijj) _ f( u_) = (B(U)Uj. 

Taking, without loss of generality, s = and linearizing about U(x), one 
obtains an equation modeling the approximate evolution of a small distur- 
bance, v. This equation has the form 

v t = Lv := {Bv x ) x - (Av) x . 

Here B(x) = B(U(x)) and A(x) is determined by the relationship Av = 
f'(U(x))v — B'(U(x))vU x . The operator L is the linearized operator about 
the wave U. 

Definition 1.1. The profile U(-) is spectrally stable if the linearized oper- 
ator L about the wave has no spectrum in the closed unstable complex half 
plane {A £ C | Re A > 0} except at A = 0. 

The next lemma allows us to narrow our search for spectrum. 

Lemma 1.1. Assuming the operator L has no essential spec- 

trum in {A G C | Re A > 0}\0. 

The lemma follows by a standard argument of [23] provided that the 
constant solutions U = U± are linearly stable. Such stability follows at once 
from the condition 

Re a(iCA(U) - \£\ 2 B(U)) < 9 > 0, 

which is equivalent to ( |++|) in the presence of f+T ) by an argument of |54| . 
Details can be found in any of ^11^31^11131 EH| Thus for the systems of our 
interest a ess (L) is confined to the left complex half plane except the origin, 
thus the only possible unstable spectrum consists of isolated eigenvalues of 
finite multiplicity. Determination of spectral stability is then reduced to 
checking that the operator L has no unstable point spectrum. The focus 
then is on the eigenvalue equation for this operator, 

Lw = Xw, 

and solutions w G L 2 with corresponding eigenvalue A such that Re A > 0. 
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The eigenvalue equation can be recast as a system of first order ODE 
(1.29) W' = A(x,X)W, W£C N . 

Because the wave U connects constant states U±, the matrix A has limits 
as x — > ±00. Thus 

A(x, A) — > A±(A) as x — > ±00. 

The idea then, is to connect the stable (resp. unstable) manifolds of the 
system W' = A(cc, A)W at ±00 to the stable (resp. unstable) subspaces 
5 ± (A)(resp. U ± (X)) of the constant-coefficient systems at each of ±00. The 
procedure we are outlining requires that the system have consistent splitting 
of the asymptotic systems. 

Definition 1.2. We say that the system kl.2ty has consistent splitting on 
SI C C if the matrices A±(A) are both hyperbolic for all A in a region SI C C, 
and there is an integer k such that the stable (resp. unstable) subspaces of 
A + and A- are both k- dimensional (resp. (N — k)- dimensional). 

The set fi is called the region of consistent splitting. In ^7] it is shown 
that the linear stability of the constant solutions is equivalent to: 

Eigenvalues /^(A) of A±(A) have nonvanishing real part 

(**) for all A with Re A > 0. 

It follows that Q, the region of consistent splitting, contains at least the 
unstable complex half plane. One important consequence of is that the 
the number of stable/unstable eigenvalues of A±(A) can be counted as A — > 
+00 along the real axis. Then given bases {4>f, ■ ■ ■ <p^} and {i>k+v ■ ■ ■ ^n} 
of the stable manifold at +00 and the unstable manifold at —00, the idea is 
to define the Evans function as 

D{X) = det(4>f, ...4>i, • • • , <^)U=o- 

A natural way to attempt such a procedure is to choose as bases for S^A) 
and U^(X) the purely exponential normal modes of W' = A±(A). Unfor- 
tunately, it is not possible to make this choice analytically with respect 
to A as some eigenvalues of A± may coalesce as A varies. The utility and 
power of the Evans function comes from the fact that this difficulty can be 
surmounted, and D(X) can be chosen to be analytic. 

The solution, due to an elegant construction of pQ, is to track volume 
forms rather than individual solutions. We associate to any collection V\, V2, ■ 
of vectors the wedge product 

Vi A V 2 A • • • A V n . 

This determines an embedding of the manifold of n-dimensional bases into 
the manifold of n-forms. More precisely this determines an embedding into 
the submanifold of n-forms expressible as a single product, the pure n-forms. 
The benefit of this approach can be seen by taking a set of n solutions 
Wi, W2, ■ ■ ■ , W n of the eigenvalue equation W' = A(x,X)W and noticing 
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that the corresponding ra-form ( = W\ A W 2 A • • • A W n solves the "lifted" 
linear ODE 

(' = A(x,\)(, 

where the operator A is determined by 

A = AWt A W 2 A • • • A W n + • • • + Wi A W 2 A • • • A AW n . 

If the collection Wi, W 2 , ■ ■ ■ , W n consists of eigenvectors of A± with cor- 
responding eigenvalues /ii,/^2, • • • ; then it is immediate that the wedge 
W\ A W 2 A • • • A W n is an eigenvector of A with corresponding simple eigen- 
value jUi + \i 2 + . . . + fj, n . Thus the volume form associated with any basis 
of <S ± (A) or U ± (\) is a simple eigenvector of A corresponding to a purely 
exponential growth or decay mode. In this lifted setting, the eigenvectors 
are simple, thus they depend on A in an analytic fashion. The construction 
demonstrates the important fact that the eigenspaces vary analytically with 
respect to a parameter even when the individual eigenvectors do not. 

For notational convenience, we follow the standard convention of associ- 
ating the full iV-volume forms with the complex numbers via the coordinate 
representation in the standard basis. That is, we write 

V 1 A---AV N = det(V 1 ,...,V N ). 

1.5.2. The Gap Lemma. The gap lemma of |17j and [22] is the key technical 
result that allows Evans function techniques for the stability analysis of 
traveling waves to be extended to the case of viscous conservation laws. This 
lemma extends the analytic framework of 1 to cases in which the essential 
spectrum of the linearized operator touches the imaginary axis, and thus 
there is no spectral gap between the essential spectrum and the unstable half 
plane {A G C | ReA > 0}. In the presence of such a gap, a standard argument 
of 8 provides a relationship between the behavior of solutions near ±oo of 
a system of asymptotically constant-coefficient eigenvalue ODEs and the 
corresponding solutions of the limiting, constant-coefficient equations. 

More precisely consider an ODE with parameter (as obtained above by 
rewriting the eigenvalue equation as a first-order system) 

(1.30) W' = A(x,X)W, 

where the differentiation is with respect to x, and A is continuous in x 
and analytic with respect to A. Moreover, suppose also that A — > A± as 
x — > ±oo. Provided that 

/-±oo 

(1.31) / |A - A±\dx < +oo, 

Jo 

then there is a one-to-one correspondence between the normal modes ^j^e^J x 
of the constant coefficient limiting system 



W' = A±(A)W, 
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where (A), (A) is an eigenvalue, eigenvector pair corresponding to A±(A), 
and solutions Wa of (jl.3Uf) having the same limiting behavior. That is 

Wp{\x) = Vfe^ x (\ + o(l)) as x -> ±00. 

The argument in [S] uses a fixed point iteration scheme depending on the 
sign of differences of the real parts of the eigenvalues /i,-. In the case of 
strict separation of the eigenvalues, a spectral gap, the fixed point is the 
uniform limit of an analytic sequence of iterates, and thus analyticity in 
A is preserved. In our case of interest there is no spectral gap; the above 
argument breaks down. The key observation in |17j was that in the absence 
of such a gap, analyticity can be preserved provided that in lieu of (|1.31j) 
the stronger hypothesis 

(1.32) |A - A ± | = OCe-"^!) as x -» ±00, 

is made. 

Theorem 1.2 (The Gap Lemma). Let A(x,X) be continuous in x and an- 
alytic in A with 

A(x, A) — > A±(A) as x — » ±00, 

at an exponential rate e~ a ' x \,a > 0, and let C — (A) an< ^ v W be analytic 
n- and (n — k) -forms associated to the complementary invariant subspaces 
of A_(A), C~ and E~ with spectral gap j3. Furthermore put tq- equal to 
the trace of A restricted to C~ . Then there exists a solution W(x,A) of the 
lifted ODE £' = A(x, \)( of the form 

W(x,\) = C(x,A)e T c-, 

where £ (and thus W) is C 1 in x and locally analytic in A. Moreover rj(x, A) 
satisfies 

(^yc(x,\) = (^) J r (A) + o( e -*M|r(A)i), 

when x < for j = 0, 1, . . . 

See IM1 EH1 El and |23| for further discussion and a proof. Appealing to 
the Gap Lemma, we thus obtain bases {4>f(x, A), . . . 4>^{x, A)} of the stable 
manifold at +00 and {4>~^ +1 {x, A), . . . 4>~^(x, A)} of the unstable manifold at 
—00. Therefore we can indeed define D(\) by 

D(A) = det(4>f, ...4%, . . . , 4>n)\x=o- 

An important feature of the construction is that -D(A) can be chosen to be 
real-valued for real A . 

Theorem 1.3. There exist bases <j>^ such that the D(X) satisfies 

D{\)=D{X). 

In particular D(\) is real valued for A € R. 
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See |63) 165) for details. A proof involves tracing through the various steps in 
the construction of the Evans function, and verifying that complex symmetry 
is preserved at each stage. 



1.5.3. The Tracking Lemma. In the calculation of the stability conditions 
(described below) it will be necessary to connect information about the 
the sign of the Evans function D(X) as A — ► oo along the real axis to 
the normalizations for the bases of stable and unstable manifolds chosen 
at A = 0. This can be accomplished by using the tracking lemma. See 

da Emms. 

Theorem 1.4 (The Tracking Lemma). For 5 sufficiently small, solutions 
w + /w~ of 

(1.33) w' = (AoOr, S) + 8(x, S))w, w G C w , 
where 5 — > is a small parameter and 

\K\ + l e l < C8, |Ao| < C. 
decaying at +00/ — 00 at rate e— x /e ax lie always within the cones 

(1.34) K_ : 



•iv 



\P s w\ < CS 



PuM 



v 



and 
(1.35) 



w 



\PuM K cs 

\Psw\ ~ 77 



respectively, for any a < liminf^^+oo a, a > limsup^.. 
there hold the following uniform growth/decay rates: 



a. Moreover 



(1.36) 
(1.37) 



\w + (x)\ 
\w+(y)\ 
\w~(x)\ 

\ w ~(y)\ 



for all x > y, and symmetrically for x < y, for any a > max x a, a < mm x a. 
(Note: Here C depends in part upon the choice of a and a. 

We note that P$ and Pu are projections onto eigenspaces corresponding to 
two different spectral groups. Rescaling in A transforms W = A(x, X)W to 
a system of the form p. 33)1 . 



1.6. Discussion. 
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1.6.1. Stability Conditions. Even though the Evans function is not typically 
evaluable, it is possible to obtain information about its zeros in the following 
way. Due to a translational eigenvalue at A = 0, D(0) = 0. One then calcu- 
lates sgnD'(O) and the sign of D(-) as A — > oo along the real axis. (Recall: 
-D(A) can be chosen to be real for real A.) Combining this information yields 
a parity for the number of unstable zeros of D(-), hence unstable eigenvalues 
for the linearized operator. We call the quantity sgnD'(0)D(+oo) the sta- 
bility index. When the signs agree, there must be an even number (possibly 
0) of real unstable eigenvalues, and when they disagree, there must be an 
odd number of such eigenvalues. Recall that complex eigenvalues occur in 
conjugate pairs, hence they do not affect the parity. Clearly then 

sgnD'(0)D(+oo) > 0, 

is necessary for spectral stability. On the other hand when the index is neg- 
ative, a positive growth rate is detected, and the wave under consideration 
is determined to be unstable. Thus the stability index is best suited as a 
predictor of instability. We remark that as the stability index only deter- 
mines the parity of unstable eigenvalues, the condition sgn D'(0)-D(+oo) > 
is not sufficient on its own to conclude spectral stability. The index yields 
only incomplete stability information; the possibilities of complex conjugate 
unstable eigenvalues and/or even numbers of unstable real eigenvalues are 
not detected by this approach. Nonetheless, the stability index serves as a 
useful starting place in stability investigation. 

1.6.2. Results. Here we describe the two main results. 

Theorem 1.5. The stability index for a strong detonation solution of Equa- 
tions M.l\) - ^I~^\j with Lax 3-shock structure has the form 

f = sgnD'(0)D(+oo) = sgn^A, 

where 7 is a constant measuring transversality of the stable/unstable mani- 
folds of the traveling wave ODE and 

(1.38) A = det(rf ,rj, [U] + q). 

Moreover for an ideal gas, the sign of the stability index is consistent with 
spectral stability in the ZND limit. 

Theorem 1.6. Strong detonations are spectrally stable for sufficiently small 
q provided the underlying gas-dynamical shock (of arbitrary strength) is sta- 
ble. 

In these theorems are right eigenvectors of the flux Jacobian, [U] is a 
vector of jumps in the gas-dynamical conserved quantities p, m, £ (density, 
momentum, total energy), and q is the vector 
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where q > represents the energy liberated during the exothermic chemical 
reaction. We note that the term A in 1)1. 38|) . which appears in the stability 
index due to the low frequency calculation of D'(0), is related the Lopatinski 
determinant itself a "stability function" for inviscid shocks. See |26j and the 
references therein. In Section 3 we detail the reduction of the equation 
A = to 

M 2 [l/p] Pe -M -1 = 0, 
the well-known instability condition of Majda |43j for inviscid shocks. Note 
that this is independent of q. We also remark that the finding T > for an 
ideal gas in the ZND limit has the implication that instability, if it occurs 
must be of "galloping" type, i.e. corresponding to the crossing of a complex 
conjugate pair of eigenvalues into the right half plane. This is consistent 
with both laboratory and numerical experiments. 

The second theorem is related to results of Liu and Ying fJJR and Li, Liu, 
and Tan j32] for various versions of the Majda model. There, full nonlinear 
stability is established for strong detonations in the Majda model when q 
is sufficiently small. In [32], the authors prove nonlinear stability for strong 
detonations in a version of the Majda model with species diffusion using 
techniques in the spirit of the Evans function. Finally we remark again, 
that by the program of [OBI ESI it is expected that spectral stability 
should be equivalent to nonlinear stability. 

1.6.3. Extensions. For the calculations in this paper, we have made the 
simplifying assumption that the equations of state are independent of the 
progress of the reaction. Though standard in the literature, this is clearly 
an idealization as the nature of the gas changes during the chemical reaction 
as pointed out in 0. One extension is to carry out the analysis in the more 
realistic setting of reaction-dependent equations of state as discussed in the 
context of the Majda model in [JJ. Also, we note that the q — > argument 
in Section ^, 7l fails since the q = gas equation is still coupled to the reaction 
equation through the equation of state. 

Another interesting direction of future study is a more detailed examina- 
tion of the effect (if any) of multiple reactants the stability index and its 
sign. In particular, the analysis of ^S] provides the geometric information 
required to evaluate (in the ZND limit) the transversality coefficient in the 
stability index for the interesting two-species reactions they consider. In 
particular, while an exothermic-exothermic two-step reaction behaves much 
as the one-step exothermic reaction we consider, an exothermic-endothermic 
two-step reaction has a richer structure |14j . 
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2. NONREACTING GAS 

In this section we consider gas dynamics as modeled by the Navier-Stokes 
equations in one space dimension; our main focus is the calculation of the 
stability index. These computations will prove useful when we shift our 
focus to detonations in the next section. The system takes the form 

(2.1) pt + (pu) x = 0, 

(2.2) (pu) t + (pu 2 + p) x = (vu x ) x , 

(2.3) £ t + {u£ + up) x = (6T X ) X + {yuu x ) x . 

The system (|2.1[) - (|2.3[1 features five unknowns (p,u,e,p,T) and three equa- 
tions. The system is completed by equations of state which incorporate the 
physical properties of the particular gas being modeled. We thus obtain a 
complete description of the fluid flow by assuming that p and T are given 
functions of p and e. The sound speed is 



c = \Jp P + P 2 PPe- 

The three differential equations (|2.1JI - (|2.3|) and the equations of state give a 
set of five equations for the five variables p, u, e, T and p. 

For some portions of the analysis, we will further assume that the gas 
under consideration is ideal and polytropic so that the specific forms of the 
equations of state are 

P{P, e) = Tpe, T(p, e) = c~ l e, 

where the constants c v and V are as in the previous section. We note that 
in this case the sound speed satisfies 

(2.4) c 2 = r e + r 2 e = (i + r)r e = 7r e , 

where j = 1 + T. 

We rewrite the system (|2.1)) - (|2.3j) in terms of the conserved quantities p, 
m, and £ , 

(2.5) Pt + m x = 0, 

(2.6) m t +i +p\ = hi — 

„t> * + (^ + p>)" = K+h=)(= 
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Rewritten once more in the form of a viscous conservation law Ut + f(U) x 
(B(U)Ux), the form of the viscosity matrix B becomes apparent 




+ d x 



0,. 



m 

m 2 1 p + p 
mp~ x {£ + p) 





- vni 





p 



6T + T e e p -^ 0T e e m + 








0T e e £ 




The Jacobian matrix A(U) of the flux / has the form 






-u 2 +p p + p e e p 
U (-E - £ +p p + p e e p ) 

To calculate the eigenvalues/eigenvectors of A, we use the device (see 
and the isentropic gas section of [HI]) of conjugating by appropriate "shift" 
matrices so that the conjugated matrix has a particularly simple form. Fol- 
lowing this procedure we obtain eigenvalues 



and right eigenvectors 



r-2 



T3 



"1 
«3 



U — C, 

u, 

u + c, 

1 

u — c 

1 

u 

2 



1 

U + C 

+ cm + - n + e 



The left eigenvectors are 



h 
h 
k 



Pp 



Pe . . Pe™ 

— e + cu + — — , • 
P 2p 

o 

+ _,-«, 1 

P 2 



Pe 

-C It, — 

P P 



Pp 



Pe 



p e u 
cu+——,c 
2p 



Pe Pe 
— U, — 

P P 
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We shall restrict our attention to a Lax 3-shock. That is, we suppose that 
the shock speed s satisfies the inequalities 

(2.8) < s < , a% < s. 

The calculations for a 1-shock follow in a similar fashion. 

2.1. Traveling Wave ODE and Linearized Equations. The traveling 
wave ODE is 

(2.9) m! = 0, 

2 \ ' / / 

m \ I I m 



(2.10) [ T+P ) = ^_ 



(2.ii) (7< £+ f>)' = (w + 



p)\P 



where without loss of generality we've taken s = 0. Each of the equations 
l|2.9|l - (|2.11jl may be integrated up once. 

(2.12) m-m_ = 0, 

m 2 \ / m 2 \ ( m 



(2 ' 13> {T +p )-{T + "^ v, 

( , 14) ty e+p) yfy e+p) )_ . («r + ,(=)(= 

The requirement for a connection that both endstates be rest points of the 
ODE, leads from (|2.12f) - (j2.14|) to the Rankine-Hugoniot conditions 

[m] = 0, 
m2 [ l /p] = -[p\, 







~P 


.P. 




.P. 



We suppose that U(x) = (p(x),fh(x),£(x)) tr is a stationary profile con- 
necting endstates U± = (p±,m±,£±) tr which satisfy the Rankine-Hugoniot 
conditions. We note that such profiles, if they exist, are transverse as an 
appropriate Wronskian is nonvanishing. Indeed global existence for connect- 
ing profiles has been shown in |2J for equations of state which satisfy the 
thermodynamic conditions of jSH]. More precisely, existence is shown when 
the equation of state is assumed to be such that the isentropes are convex 
and do not cross in the pressure-volume plane. We note that a polytropic, 
ideal gas satisfies this condition. 

Then linearizing about this profile, we find equations for the evolution of 
small perturbations (p, m, £). These equations can be written in the general 
form 

w t + (Aw) x = {Bw x ) x , 



2-1 
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where w = (p, m, £) tr and the matrices A and B depend only on x. More 
precisely we find 

(2.15) p t + m x = 0, 

(2.16) m t + (a 2i (x)p + a 2 2(x)m + a 2 3(x)£) x = (b 2 i(x)p x + b 22 (x)m x ) x , 



(2.17) £ t + (a 3 i(x)p + a 32 {x)m + a 33 (x)£) x = 

{hx{x)p x + b 32 (x)m x + b 33 (x)£ x ) x , 

where the coefficient functions ctij(x) can be expressed in terms of the entries 
of the flux Jacobian and derivatives of the entries of the viscosity matrix. 
The terms bij{x) correspond to the ij-entries of the viscosity matrix B. 
In both dependence arises from evaluation along the known profile 

U{x). 

As our interest is in spectral stability of the profile, we focus on the 
eigenvalue equations corresponding to (|2.15|) - (|2.17|) . They are 

(2.18) Xp + m' = 0, 

(2.19) Am + (a 2 i(x)p + a 22 {x)m + a 23 (x)£)' = (b 2 i(x)p' + b 22 {x)m')\ 



(2.20) \£ + (a 31 (x)p + a 32 (x)m + a 33 {x)£)' = 

(b 3 i(x)p + b 32 (x)m + b 33 (x)£' 
and the corresponding limiting system as x —* ±oo takes the form 

\p + m! = 0, 
Am + (af t p + af 2 m + af 3 £)' = {b^p + bf 2 m')' , 
X£ + (af lP + af 2 m + af 3 £)' = {bf lP ' + bf 2 m' + bf 3 £')'. 
We make the invertible change of variables as in 



V 






H 










1 


o > 






&21 


b 22 





II 


m 1 


b 3 i 


b 32 


^33 J 




V s ) 



Thus 



and 



C 

m, 

b 2 iP + b 22 m, 

b 3 \p + b 32 m + b 33 £, 



-b 2l %2 b£ 



c~ L - I 1 

b 2i b 3 i b 31 b 22 - 633 1 b 32 - b 3 £ b 2 l 631 633 1 
We can rewrite the eigenvalue equation in z-coordinates as 

{B(C- x zj)' = (AC-h)' + AC -1 *, 
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or more explicitly as 

(2.21) = z[ + X( y -b^b 22 z 1 + b^z 2 ), 

(2.22) z 2 = (PiZ! + p 2 z 2 + p 3 z 3 )' + Xz u 

(2.23) z 3 = (rjizi + r] 2 z 2 + r] 3 z 3 y + Xg(z x , z 2 , z 3 ), 
where the linear function g is given by 

(2.24) g(z 1 ,z 2 , z 3 ) = [b 2 lb^b 3 ib 22 - b^b 32 )zi - (o^o^&si)^ + fc£z 3 , 

and the (3 and n coefficients can be calculated in terms of the and bij . The 
exact form these coefficients is not used below, so we omit the calculation. 
From (|2.21|) - (|2.23j) it is a simple matter to recast the eigenvalue equation as 
a first order system of the form 

Z' = A(x, X)Z, Z = ( Zl ,z 2 , z 3 , 4, z' 3 ) tr , 

with a corresponding limiting system 

Z' = A±(A)Z 

at each of ±oo. The matrix A takes the form 



A(x, X) 



t 


-Xb^lb 22 


Xb 2 l 








\ 













1 



















1 




A + # - Xb^b 22 


(3' 2 + Xb 2 l 


ft 

V 3 + 




ft 


\ 


A 5 i 


A 52 


m 


V3 J 



where 
and 



A51 =Vi- ^ b 2 i b 22 + A(6 21 1 6 3 3 1 63i622 - 633^. 



32 J 



A 5 2 = TJ2 + Xb 21 ~ ^ b 33 b 21 b 31- 

We verify the consistent splitting hypotheses, without loss of generality, in 
the original ^-coordinates. The characteristic equation has the form 

(XI + nA± - n 2 B±)v = 0. 

We obtain then, a sequence of lemmas. 

Lemma 2.1. For Re A > the matrix A±(A) has eigenvalues 

M ±(A),/4(A) < < M3 fc (A),//±(A), At ±(A), 

funi/t ordering referring to real parts). The eigenspaces S (A) and ^/^(A) 
associated with the eigenvalues nf (A),/x 2 (A) and // 3 (A), /i 4 (A), // 5 (A) re- 
spectively depend analytically on X 

Proof. As established in the introduction, the number of positive/negative 
roots is constant for Re A > 0, so that roots can be counted as A —* +00 
along the real axis. There is one root with fi ~ A and four roots with 
[i ~ A 1 / 2 . See the Appendix for more details. □ 

Moreover in a neighborhood of A = 0, a bifurcation analysis yields: 
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Lemma 2.2. For each j, there are analytic extensions of (A) to a neigh- 
borhood N of X = 0. Moreover there are analytic choices of individual 
eigenvectors Vj (A) corresponding to fij (A) m TV. 

Proof. When A = the characteristic equation reduces to {fiA± — fJ?B±)v = 
which has a triple root at zero. Nonzero roots must satisfy € a(A~ l B) 
which has two nonzero eigenvalues. One of them switches signs at ±oo. The 
zero roots bifurcate analytically from zero, for linearizing about (A, fi) = 
(0, 0), we obtain 

(\I + fj,A±)v = 0, 

we find that on the — oo side 

^3,4(A) = -^+0(A 2 ), 

a l,2,3 

and 

V 2,3A = r 1,2,3+ °( A )- 

While on the +oo side 

^4, 5 (A) = -^+0(A 2 

and 



a l,2,3 



Then V- correspond to v ■ when the limiting eigenvalue equation is written 



^3,4,5 = r 1,2,3+ 

r± 
j 

as a first order system. □ 
Lemma 2.3. There are choices of bases 

Bf(A) = {0±(A),0±(A)} 

and 

B±(A) = {^(A),0±(A),0±(A)} 

ofS :il (X) andU zil (X) which are analytic with respect to X in iVU{ReA > 0}. 
In the neighborhood N , they satisfy 

B + (X) = k + (X)V+(X)AV + (X), 

and 

B U {X) = MA)Kf(A) A V 4 {X) A Vf(X) 

where Vp are as in the previous lemma and k±(X) are scalar functions such 
that k±(0) = 1. 

Proof. The proof follows from the previous lemma and a standard (nontriv- 
ial) result of matrix perturbation theory j^H] . □ 

Finally using the gap lemma (Theorem ll.2|) . we obtain 
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Lemma 2.4. There are bases Bs(x,X) and By(x, A) of the spaces of solu- 
tions of the eiegnvalue equations decaying at x = ±00 which are tangent to 
S + (X) as x — > +00 andU~(\) as x — > —00. That is 



and 



<S + (A) = lim span jBs(x, A), 



U (A) = lim span£>^(x, A). 

X— > — 00 



A word on notation is in order. Working again in z-coordinates, we refer to 
the elements of the bases Bs(x, A) and Bu(x, A) by Zf{x, A), so that 



Z t ~ ( z tj ' z 2,j ' z 3,j > z 2j > z 3j' )* r > 



and we denote by with a solitary subscript the first three components of 



Zf. Thus 



Z± = (zf .,zf .,zf .) 
J ^ lj' 2,J : 3,3 J 



2.2. The Evans Function. 
Definition 2.1. XYie Evans function is 

£>(A) = det(Z+, Z+, Z 3 -, Z 4 -, Z 5 -)| x=0 = det 





Z l,2 


z l,3 


Z l,A 


z iA 


Z 2.\ 


Z 2,2 


z 2,3 


Z 2,A 


z 2,5 


z 3,l 


Z t_2 


z 3,3 


Z 3,A 
/ 


z 3,5 


z 2,l 


Z 2,2 


z 2,3 


z 2,4 


z 2,5 


Vt'i 


z 3,2 


_/ 
z 3,3 


_/ 
Z 3A 


Z 3,'J 



\x=0 



As usual we are free to put at A = 



CU X , 



(2.25) zj = z, 
and at A = we also set 

(2.26) 4(+°°) = °> z 3(-°°) = Crf, z 4 (-00) = Cr 2- 
2.2.1. Calculation of D'(0). 

Proposition 2.5. T/ie Evans function D(\) satisfies D(0) = and 
sgnD'(O) = sgn det (rr,r^, [[/]), 

w/iere 

/ Z 2T Z 2^2 

7ivs =[-■(■ 7 + 
V ^3,1 ^3,2 

which measures transversality of the intersection of stable/unstable mani- 
folds in the traveling wave ODE. 
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Proof. That -D(O) = follows immediately from (|2.25j) . Using the Leibniz 
rule to compute -D'(O), we find 

D'(0) = det(d x Z+, Z+, Z 3 -, Z 4 -, Zs)\ x=0 + ■■■ 

+ det(Zf , Z 2 , Z 3 , Z 4 , (?a-^5 )U=o- 
We combine the two nonzero determinants above 

(2.27) D'(0) = det(Z+,Z+,Z 3 -,Z 4 -,Z)| x=0 , 
where 

z = 5 A (z 5 - - Z+). 

Differentiating the eigenvalue equation with respect to A leads to the equa- 
tions satisfied by z. 

(2.28) B{C- l z)' = AC- 1 z + [U]. 

Also at A = 0, the eigenvalue equations simplify considerably for j = 1,2,3,4 
to (omitting ±) 

o = 4,p 

z 2,j = (Pl z l,j + & Z %3 + 03Z3,j)', 
Z 3,j = (VlZlJ 1 + V2 z 2,j + V3 z 3,j)', 

which can be integrated up using the boundary conditions supplied by (|2.25|) 
and (t2~26l) . Thus 

(2.29) B(C- l z+)' = AC- l z+ j = 1,2, 

(2.30) B{C- l z^)' = AC- 1 z 3 -a^r^, 

(2.31) B{C~ l zl)' = AC" 1 zl - a 2 r 2 . 

The first equation of each of (|2.28j) - (|2.31|) allows a simplification in the first 
row of the determinant (|2.27|) . Namely, 

z i = ~[p], 
zfj = j = 1,2, 

z i,s = «r( r r)i' 

Z ~ A = a 2 {r^)i. 

The second equation of each of (j2.28j) - (|2.31|) allows a row operation to sim- 
plify the fourth row, 

z 2 = 0izi + (3 2 Z2 + /?3^3 + [m], 
z 2,'j = 0izt,j + foztj +fcztj, 3 = 1, 2, 

z 2,3 = @l Z l,3 + 02*2$ +^3^3 - ar( r D 2, 

z 2A = fcz^ + fcz^ + fcz^-ciz (r^) 2 , 
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while the third equation indicates that a row operation will simplify the fifth 
row of (l2~2 



z 3 



"3,3 



"3,4 



1,2, 



rjizt + ij 2 z 2 + mh + [£], 



where we've used {r- )j to denote the zth component of r ■ . Putting these 
operations together, (|2.27j) simplifies to 



(2.32) 



det 



/ 





V o 



"1 








- a i_( r i_h 

~ a l( r l)3 



«2 

-a 2 (r 2 h 

-a)2-(r 2 ) 3 



-IP] \ 
* 

[m] 



From (|2~32j) it follows that 
D'(0) 



a 1 a 2 det 



21 2 2 I det(rf ,r 2 





: 



















Z 2,l Z 2,* 

Z + Z + I 
z 3,l z 3,2 J 

From the shock inequalities (|2.8j) it follows that a^a^ > 0, and the propo- 
sition is proved. □ 

2.2.2. Large A Behavior. We appeal to the Appendix to determine sgnD(A) 
as A — > +00 along the real axis and complete the calculation of the stability 
index. In the Appendix the calculation is carried out for abstract "real 
viscosity" systems of the form: 

U t + F(U) X = (B(U)U X ) X , 

where 


and 

n,/GR n " r , v,geR r , &i € R rx(n " r) , b 2 £ R rxr . 
Noting the abuse of notation, for this section we adopt the notation of 
Appendix. Then from Lemma lA.21 it follows that for A real and sufficiently 
large, 

D(X) + 0. 

Furthermore we can relate sgn.D(A) for A large to our normalizations at 
A = by Lemma lA.31 We find for A real and sufficiently large that 

(2.33) sgn£>(A) = sgnS+ det(vrW + , eS + ) det(7rW~)| A=0 , 

where D is the Evans function computed in the original w-coordinates. D 
and D differ by a nonvanishing real factor. In (j2.33fl . §> + is a basis for the 
one-dimensional stable subspace of A = An — A^b^bi, where Aij are the 
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entries of the Jacobian matrix of F, in particular (An, A12) = df(U). In this 
case A is simply the particle velocity (the "original" u) and is thus negative 
for a 3-shock. W 1 * 1 denote bases for decaying solutions at each of ±00, and 
7T, e denote projection and extension respectively. 

Working now in z-coordinates and following the discussion of the Appen- 
dix, the form of (|2,33|) simplifies considerably, and D(\) satisfies 

Lemma 2.6. For A real and sufficiently large 

(2.34) sgnD(A) = sgn 7w5 (S + ) 2 det(rf , rj, U x ). 

2.3. The Stability Index. Therefore combining Lemma 12.61 and Proposi- 
tion [23] we have 

Proposition 2.7. The stability index for a viscous Lax 3-shock is 

sgnZ?'(0)D(+oo) = sgndet(rf,7YT, [U]) det^f,^, Uxl-oo), 

where [U] = ([p], [m], [£]) tr is a vector of jumps and U x = (p x , rh x , £ x ) tr ■ 

We remark that the stability index is unaffected by change in the viscos- 
ity matrix. That is, the index agrees with that of an artificially parabolic 
system at least in the presence of (l+l)- (l++'+l ). Further we note that the cal- 
culation of D'(0) captures low-frequency information. In this setting, this 
corresponds to "inviscid" behavior or the stability of shocks solutions of the 
Euler Equations, 

Pt + (pu)x = 0, 

(pu) t + (pu 2 + p) x = 0, 

£ t + u{£+p) x = 0. 

and thus A takes the form of a Lopatinski determinant, familiar from the 
stability analysis of such shocks. See [2^1 f° r the calculation of this deter- 
minant for the multi-dimensional Euler equations. In the weak-shock limit 
we note that 

[U] ~ U x ~ r^, 

which implies that the stability index satisfies 

sgnD / (0)D(+oo) ~ sgndet(r^, r^, r^) 2 = +1 

in the weak shock limit consistent with stability. Moreover in the ideal gas 
case, combining the nonvanishing of A which is well known for an ideal 
gas, for example see [33] or [221, the global existence result of |21| which 
guarantees transversality of connections, and nonvanishing of D(+oo), we 
can then conclude consistency with stability for shocks of any strength in 
the ideal gas case. 

We remark that in the more general case, some rudimentary knowledge 
about the connecting orbit is neccessary to evaluate 

sgn det{r^,r^,U x ), 
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namely the direction of U x \- <x . 

2.4. Isentropic Gas Dynamics. We note that the in case of isentropic 
gas dynamics, 

(2.35) p t + (pu) x = 0, 

(2.36) (pu) t + (pu 2 + p) x = (vu x ) x , 
where the pressure satisfies 

p = p(p), c 2 =p\p) > 0, 

an analogous but simpler calculation leads to the stability index 

f = sgnZ/(0)£>(+oo) = sgn7 i 2 det(S+) 2 det(rf , [U]) det(r~ ,U X ) 

= sgndet(rf, [U])det(r~ ,U X ), 

for a 2-shock. Here [U] = ([p], [m]) tr is a vector of jumps in the gas- 
gynamical conserved quantities, U x = (p x ,m x ) tr , r{" = (1,U- — c_)* r is 
the outgoing right eigenvector, and the term 7, is a transversality coeffi- 
cient. In this case, evaluation of the index is straightforward since the phase 
space for the traveling wave ODE is one-dimensional. The phase portrait is 
shown in Figure l2~Tl Then sgndet(rj~, [U]) det(rj~, U x ) becomes simply 

« I » p 

p+ p- 

Figure 2.1. The Phase Portrait 

sgn[p]p x (u~ — cJ) 2 = +1, consistent with stability, 

since both [p] and p x are negative. Moreover, we note that for one-dimensional 
isentropic gas dynamics, Humpherys j2E] has shown that there are no un- 
stable real eigenvalues for shocks of arbitrary strength. 
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3. Reacting Flow 

In this section we extend the analysis of the previous section to our main 
interest: strong detonation waves, which are particular traveling wave solu- 
tions of equations (|l.l)) - p.4j) . We assume "ignition temperature kinetics," 
that is, (p vanishes for temperatures below some ignition temperature and is 
identically 1 for some larger value of T. We also assume that the pressure 
p and the temperature T are given functions of p and e, and thus are inde- 
pendent of the progress of the reaction. At some points in the analysis we 
further specify that the gas is ideal and polytropic, so 

p = Tpe, T = c~ l e. 

Here e is related to e, the specific internal energy not due to reaction, by 

e = e + qY. 

Using this relation, the third and fourth equations can be combined, and 
the energy balance equation can be rewritten as 

P (y + C )) + ( K PU ^~2 + e ) + np ) = (^)* + qkpY(p(T) + (uuu x ) x . 

Using this form of the energy balance equation, we note that the system fits 
in the general framework (|3,87|) - (|3.88|) discussed below with U = (p, m, S) tr , 
z = p Y, Q = q = (0,0,q) tr , $(U) = kip(T), D 1 ^) = d, and D 2 (U,z) = 
(—^,0,0). We remark at this point that since the flux / in the kinematic 
variables U is as in the Navier-Stokes case, the left and right eigenvectors 
lj,rj of the flux Jacobian A(U) and their corresponding eigenvalues aj are 
precisely as calculated in the previous section 

In accordance with our previous analysis, we assume that the strong det- 
onation has Lax 3-shock structure, so 

al^ < s < a J , s < a% . 

We will without loss of generality set s = 0, so the reaction front is station- 
ary. The shock inequalities imply that u± < 0, so that fluid particles cross 
the reaction front from right to left, or alternatively the front is "moving" 
to the right connecting an unburned state at +oo to a completely burned 
state at — oo. Without loss of generality we normalize p+ = 1 so that the 
total reactant variable z satisfies 

(3.1) z+ = 1, z_ = 0. 

We also assume that the endstates are such that the temperature on the 
unburned side is below ignition so ip = and that the temperature on the 
burned side is sufficiently large so that <p = 1. We'll see that this first 
assumption is necessary to ensure that the end state at +oo is a rest point 
of the traveling wave ODE. Thus ip satisfies 

(3.2) <p + = 0, ip. = 1. 
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3.1. Traveling Wave ODE and Linearized Equations. The (s = 0) 

traveling wave ODE is 

(3.3) m a 



(3.4) 
(3.5) 
(3.6) 



= 0, 

to 2 

Hp 

P 



p 

( mz 
V P 



P 



m\ I to 

11 '"7 (7 



(dz x ) x + 



P 



-Pa 



k(p(T)z. 



Fixing a state at — 00, the first three equations above can be integrated once 
to 



to — to_ 

2 

to 



0, 



P 



m , ~ 



\-P Hp 



m 
v I — 



P 

0T X + 



qpd ( - 

VP 



+ — — 



p 



m 



P 



Requiring the state at +00 to be a rest point of the ODE yields from the 
first two equations above the familiar Rankine-Hugoniot conditions 



(RH1) 
(RH2) 



TO 



TO 



P 



+ P 



0, 
0. 



where the brackets as usual indicate the difference between the state at +00 
and that at —00. In the case of the third equation, we find 



m£ mp 
P P 



m£ mp 
P P 



mz 



since £ = £ + qz. Moreover from (|3.1|) . the third jump condition becomes 



(RH3) 







m£ mp 
P P 



qu + . 



Linearizing about a profile (p,m,£,z), we find 

(3.7) Pt + m x = 0, 

(3.8) m t + (a 2 iP H h a 2 3£) x = (b 2 ip x + b 2 2m x ) x , 

(3.9) £t + (&3ip-\ \-a 33 £) x = (b 3 ip x -\ \-b 33 £ x ) x + qkl, 

(3.10) zt + (vip + v 2 m + v±z) x = (dz x ) x + (dp x ) x - kl, 
where and \j are as in the previous section, and 

I = {l z z + l p p + l m m + le£}, 
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with 



and 



lz(x) : 
1 P (X) : 

k{x) ■■ 

vi{x) 
v 2 (x) 
f 4 (x) 



<p'(f)(f p + f e e p ), 
P 



-uY 



dzp x 



z 

~ 3 
P 



dp x 
P 



We note that due to the structure of <~p and 1)3.2(1 . it follows that l p , l m , and 
lg vanish at both ±00 while l z+ = and l 2 _ = 1. The equations (|3.7|) - (|3.1U|) 
can also be written in the more compact form 

(Bw'Y = (Aw)' + w t + qkg(w,z), 

(dz')' + (dw'Y = (V w wY + {V z z)' + z t -kg(w,z). 

In an abuse of notation we've written d to stand for both the 1x3 matrix 
(—^,0,0) and the (1, l)-entry of that matrix. The meaning will be clear 
from the context. 

Under the invertible change of coordinates 

/ Ci \ / 1 \ 

C2 



(3.11) 



V c 4 / 



hi 
hi 
V d 



h2 
h2 






&33 








d ) 



( p \ 

m 
£ 

V * J 



c 



the eigenvalue equation corresponding to 1)3.7(1 - 1)3.10(1 takes the simple form 



(3.12) 





= C! + A(- 


-b 21 1 h 2 ( 1 +b 2 - 1 1 ( 2 ), 


(3.13) 


C2 


= (ACi + 


••+/?3C 3 )' + ACi, 


(3.14) 


& 


= (mCi + 


• • + mCs) 1 + Xg(0 + gki-( 


(3.15) 


It 

G 


= + 


■ ■ + 9aCaY + A/t(C) -kl-(. 



We have used the notation £ = (£1, ( 2 , £3, £4) and C = (Ci , C2 , (3)- The linear 
functions g and h take the form 



and 



9(0 = (6^6^631622 - 6i 3 1 632)Ci - (6i 3 1 6 21 1 6 3 i)C2 + 633^3, 
fc(C) = d-%lb 22 dd - d-%ldC 2 + cr 1 ^. 
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The coefficients (3j and rjj are as in the previous chapter, while the 6j depend 
on d,d,bij and Vj. Also I ■ ( = + ■ ■ ■ + Z4C4) where lj depends on 

lp, lm, hi lz,bij,d and d. 

The change of coordinates matrix C has block structure respecting the 
division of the variables into gas dynamical variables w = (p,m,£), and 
reaction variable z. Thus we have 



c 



Cms 



d 







d 



w 

z 



where C^g denotes the change of variables in the nonreacting gas dynamics 
case considered in the previous section. The inverse of the matrix C also 
respects this structure and is given by 



C 



-1 



d 1 b 2 lb 22 d 



n- 1 



-d 



-%ld 







dr 



We note that when A = 0, by a substitution of (|3.15|) into (|3.14|) through the 
term I ■ £ which appears in both equations, we recover a third equation in 
which every term is differentiated; this means that we can integrate up the 
first three equations subject to appropriate boundary conditions. This key 
fact is what will allow for a simplification of the Evans function determinant 
via row operations in the calculation of D'(0). The eigenvalue equations 
(|3.12|) - (|3.15j) . can be written as a first order system of the form 



/ \tr 
4) ) 



(3.16) Z' = A(x,\)Z, Z = (Ci,<2,C3,C4,Ca,C3»C 
with corresponding limiting constant-coefficient system 

(3.17) Z' = A±(\)Z. 

The characteristic equation in original (w, z)-coordinates (with w = (wi, W2, ^3) 
is 



tr \ 



M± 



w 

z 



where the characteristic matrix M± has block structure 



with blocks 



M. 



ui 



V M031 





/ Mf 








\ M l l 


M l l ) ' 




A 









- M 2 ^2i 


A + pa 2 2 - f-i b 2 2 


pa 2 3 


- M 2 ^3i 


pa 32 


~ M &32 


A + /i«33 - 



Ml r 




36 



GREGORY LYNG AND KEVIN ZUMBRUN 



and 



MS 



M 



(fxvi - fi 2 d /iV2 0) , 

± — A — Li 2 d + /xv 4 + /c(^±. 

As usual we take advantage of consistent splitting and count the number 
of stable/unstable roots in the limit A — > +00. Due to the incomplete 
parabolicity of the viscosity matrix B, we expect one root which scales as 

fx ~ flX, fx ~ 1, 

hence the matrix M± takes the form 

f\ flX 



fxXa^i 
jxXa 3 i 



fx 2 X 2 b 
A 2 A 2 &3i 



21 



A + /iAa 2 2 - fl 2 X 2 b 2 2 



fxXa 32 - fx z X z b, 



:>>> 





flXa.23 
X + fla 33 - fl 2 X 2 b 3 3 






-qktp± 



\ 



fi 2 X 2 d 



flXY± 







iV / 



\ /iAu ± i± 

where 

N = X — fl 2 X 2 d + flXu± + k<p±. 

Whereupon dividing each row by the highest power of A, we find that in the 
limit A — > +00, the roots satisfy a block triangular system 



(3.18) 



/ 


1 


f< 





\ 




f \ 




b2i 


b22 








U>2 







631 


&32 


^33 





W 3 





V 


d 








dj 


\*J 


w 



where the upper left-hand block is easily recognized as the gas dynamics 
block which appears in the Appendix. Thus there is one root 

-fT 1 G <t(u±), 

and therefore since we are considering detonations with 3-shock structure, 
there is one unstable root at each of ±00. The remaining roots scale as 

£A 1/2 , 



fx 



1. 



Thus, upon substituting and dividing as before, we obtain a different block 
triangular system 



(3.19) 











° \ 








-fl 2 b2i 


-fx 2 b 2 2 + 1 








W2 







fl 2 b 3 i 


-ft 2 b 3 2 


-fi 2 b 33 + 1 





W 3 







-fx 2 d 








-fx 2 d+l J 






w 



It follows immediately from (J3.19|) that w\ 



b22 









b 3 2 


&33 













\d) 





0, and thus roots must satisfy 

M 
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that is, fl 2 must be an eigenvalue of the block diagonal matrix 





/ b 22 








(3.20) 


b 3 2 


&33 







V o 





d 



This yields 3 stable and 3 unstable roots at each of ±oo. We note here that 
due to the block diagonal structure of (j3,20|) . there is a 2-1 split at each 
infinity of kinematic versus reactive roots. Summarizing, we've found that 
there are 4 stable and 3 unstable roots at each of ±oo as long as Re A > 0. 

Consistent splitting breaks down at A = 0. We look first at the x = +oo 
case. In this case due to (|3.2[) . the characteristic equation has lower block 
triangular structure when A = 0. Thus the characteristic matrix has the 
form 










A* 










- A« 2 &2i 




- I^b 22 


Ai0 2 3 







- Ai 




- Ai 2 &32 


/"a 33 - fi 2 b 33 





fjtu+ 


-^d 









— fi 2 d + fj,u+ 



\ 



J 



From our analysis of the nonreacting case, we know that the stable roots 
corresponding to the upper left-hand kinematic block do not vanish at A = 0, 
thus they correspond to fast kinematic modes. On the other hand, looking 
at the reaction block, 

—fi 2 d + = 0, 

it's clear that there are solutions \x = corresponding to a slow unstable re- 
active mode and fi = a fast stable reactive mode. A bifurcation analysis 
as in the previous section shows that the roots have analytic extensions in 
a neighborhood of A=0. 

In the x = — oo case, we find when A = that, due to (|3.1|) . the charac- 
teristic matrix is upper block triangular 

/ n 

/ia 22 - /i 2 ^22 /ua 2 3 

fia 32 - n 2 b 32 /ia 3 3 

\ I -\i l d + /j,u + k ) 

From our previous analysis of the nonreacting case, we know that two of the 
unstable roots from the gas block vanish at A = 0. By block structure, the 
corresponding vectors have the expansion 





fj,a 2 i - fJ- 2 b 21 
fj,a 3 i - fi 2 b 31 



A3 






-qk 



\ 



+ 0(A), J = 1,2, 



where r- is an eigenvector of the flux Jacobian A. The roots coming from 
the reaction block must satisfy 



U- 

~2d 



u 2 . + idk 



2d 



7^ since d, k > 0. 
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This implies that all reactive modes are fast modes on the x = — oo side. 
We note here that this discussion shows that all the reactive modes of our 
interest (stable at +oo/unstable at — oo) are fast modes, and thus asymp- 
totically vanish in both kinematic and reaction components. We note here 
that the inclusion of the differing block structure at each of ±00 corrects an 
omission in the abstract development considered in |65j . 
Finally, applying the gap lemma we obtain bases 

{Z+,Z+,Z+} and {Z 4 -,Z 5 -,Z 6 -,Zf} 

spanning the stable manifold at +00 and the unstable manifold at —00. We 
use the notation 



Z f 



(Ci, j ) ctj > (Tj ' Cy > (Tj > c^j j cTj ) tr ' 



±' /-±' /-±'\tr 



when we need to indicate the components. Also we note here that we will 
follow our standard convention and identify /Zj with the derivative of 
the profile since it satisfies the linearized system and decays at both ±00. 
We'll also let Z^ and Z^ correspond to fast kinematic and reactive modes 
respectively. On the —00 side we'll let Z± and Z^ correspond to slow 
kinematic modes while Zq is a fast kinematic mode. Note that implicit in 
this assignment is the assumption that the profile approaches the burned 
endstate parallel to the reactive mode. This assumption is generic in the 
case where the reaction is slower and actually occurs in the ZND limit as 
evidenced in our discussion of existence in Section 1. 

3.2. The Evans Function. 
Definition 3.1. The Evans function is 
(3.21) D(X) = det(Z+ Z+, Z+ , Z±, Z~ , Z 6 , Z?)\ x=0 . 

We make the usual normalizations at A = 0. Namely we put 

\ 

, + 

(3.22) 



^3,7 

V C77 



c 



( Px 



\ *» J 



and the other fast modes satisfy 

Z. 

while the slow modes satisfy 



(3.23) 



2 (+°°) 

6 ("OO) 



Z 3 +(+oo) = 0, 
0, 



(3.24) 



/ Cm \ 

C2,4 
C3,4 

V C4-4 / 



c 



( Ci> \ 

^2,5 
^3,5 

V C4-5 / 



c 
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3.2.1. Calculation of D'(0). 

Proposition 3.1. The Evans function D(X) satisfies D(0) = and D'(0) = 
7dA where 

(3.25) A = det(rr,r 2 -, [[/]+<?), 
and 

( C^l ^2,2 (£,3 C2 ^ 

(3.26) ld = det 



C3J v',.2 Ci.3 V;.<; 
C44 Ca2 S4,3 C4 



s4,l 

\ C4,l C4,2 £4,3 Q,6 / 

measures transversality of the stable and unstable manifolds in the traveling 
wave ODE. 

Before beginning the proof we remark that in contrast to the case con- 
sidered earlier 7^ is not extreme, that is, it involves fast modes from both 
infinities. 



Proof. As usual D(0) = follows immediately from the normalization (|3.22|) 

chosen for the basis elements. Applying the Leibniz rule, 

(3.27) 

D'(0) = det(8 x Z+, Z+, Zf)| x=0 + • • • + det(Z+ . . . , Z 6 " d x Z^)\ x=0 , 

and combining the two nonzero determinants in the above equation, we 
obtain 

(3.28) D'(0) = det(Z+,Z+, Z+ , Z 4 " , Z 5 ~ , Z % , Z)\ x=0 , 

where Z = Z~ — Z + = d\Zj — d\Z±. For clarity we perform the neces- 
sary manipulations in the original (w, z)-coordinates and then translate the 
results to the ("-coordinates. Thus we write down the general form of the 
eigenvalue equation 

(3.29) (Bw')' = (Aw)' + Xw + qkipz 

(3.30) (dz')' + (dw')' = (V w w)' + (V z z)' + Xz - kcpz. 

When A = 0, we can make a substitution from ()3.3(J|) into l)3.29j) so that the 
first equation has the form 

(3.31) (Bw'Y = (Aw)' + q(-(dz')' - (dw')' + (V w w)' + (V z z)'), 

in which every term is differentiated and hence (|3.31j) can be integrated to 

Bw ' —Bw± = Aw—Aw±+q(—dz'+dz±—dw'+dw±+V w w—V w w±+V z z—V z z±) 

with the db subscripts indicating boundary conditions to be supplied by the 
normalizations (|3.22)) - ([3.24[) . Therefore the fast modes satisfy 

(3.32) Bw' + qdz + qdw' = (A + qV w )w + qV z z, 
while the slow modes satisfy 

(3.33) Bw' + qdz' + qdw' = (A + qV w )w + qV z z — ajrj . 
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On the other hand (w^, z^) satisfy the variational equations at A = 

(3.34) (Bw^)' = (Aw*)' + U x + qki P z ± 

(3.35) (dz ± ')' + (dw+'y = (V^y + (V^)' + z x - ktpz* . 

We make the same substitution from (|3.35fl into (|3.34j) so that every term 
is a derivative. Then we integrate (w + , z + ) from x to +oo, 

(3.36) Bw +> + qdz + ' + qdw + ' = (A + qV w )w + + qV z z + + U -U + +z- z+, 
and (w~,z~) from — oo to x, 

(3.37) Bw~' + qdzT* + qdw~' = (A + qV w )w~ + qV z z~ +U-U- + Z-Z— 
It follows by subtracting (|3.36|) from (|3.37l) that (w, z) satisfy 

(3.38) Bw + qdz + qdw = (A + qV w )w + qV z z + [U] + q 

Translating this information to (-coordinates, we have for fast modes (kine- 
matic and reactive) j = 1,2,3,6 

o = Cu 

Note that the equation for £4 ■ is unchanged. We don't get any simplification 
from that equation. Also when j = 4, 

= Cl,4-ar( r l~)i 

c^4 = PiCia h — i- /%C3,4 - ar( r r)2 

c'3,4 = 771(1,4 h — 1- r? 3 c3,4 + qK-ci^ + ^iCi,4 h — h&iCm) _a r( r r)3) 

and similarly for j = 5, 

= Ci,5- aJfcDi 

C 2 ,5 = ACl,5 H h /?3C3,5 - a 2 ( r 2^)2 

C3,5 = mCi^H h r/ 3 C 3 ,5 + ^(-(1,5 + 01&.5 H h 6*4(4,5) - 0^(^)3. 

Finally ( satisfies 

= Ci + [p] 

I2 = PiCi + • • • + /%& + H 

('3 = V1C1 + ■ ■ ■ + mCs + qH-(' 4 + 9iCi + • • • + 4 a) + [£} + q 
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These equations allow us to perform row operations simplifying rows 1, 5, 
and 6. Then rearranging rows and columns we find that 



£>'(0) = det 













\ 












<h r l a 2 r 2 i U ] + 9 
















(2,2 


^2,3 


(2,6 


* # * 


<& 


(3,2 


^3,3 


£3,6 






(4,2 


C4,3 


(4,6 




«1 


c +' 

S4,2 


c+' 

M,3 


^4,6 


* * * / 



from which the result follows. 



□ 



3.2.2. Large A Behavior. To finish the calculation of the stability index, it 
remains to determine 

sgnD(A) as A —> +00, real. 

Here block triangular structure plays a key role. In particular it allows for 
certain "diagonal" kinematic and reaction components of eigenvectors to 
be separately analytically specifiable. This means that we will be able to 
calculate the sign in two pieces. For the gas dynamical piece the analysis of 
Appendix C applies, while the reaction piece can be treated directly. 

For ease of comparison and consistency with our previous analyses, we 
shift now to the notation of Appendix C. Thus, (u, v) represent the gas 
dynamical variables, A = An — A^b^ b\ where Ay is the tj-entry of the 
flux Jacobian and b\, 62 are the nonzero blocks of the matrix B. Recall that 
in this case A is simply the particle velocity (it in the usual notation). For 
3-shocks, this quantity is uniformly negative. We also use z to represent the 
reaction variable. 

Proposition 3.2. For real A sufficiently large, 

(3.39) sgn£>(A) = sgndet(S+) det(yrW + , e§+) det(yrW") x z+zf / 

kinematic only j'=l,2,4,5,6 reaction 

where ir denotes projection of the kinematic variable W = (it, v, v') onto 
(it, v) components, and S(x) is a real basis of the stable subspace of A, and 
EU := (it, — 1> 2 



-b n 1 b\u) denotes extension. 



Proof. For real A sufficiently large, from the block triangular structure in 
H3.18j) and (|3.19|) . it follows that at any (fixed) x the stable/unstable mani- 
folds of the frozen eigenvalue equation are spanned by vectors of the forms 



(3.40) 



/ u \ 

-b^biu 



( \ 

v 

Tft\ 1/2 v 

* 
* 



(kinematic) 
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and 



(3.41) 








1 



(reaction) . 



Therefore the Evans function satisfies 

r / o o 



D(A) - det < 



Vl 

-fl 1 X 1 / 2 v 1 
* 
* 



V 2 

-jl 2 X l / 2 v 2 

* 



u 

-b^biu 

* 

* 



where and are as in the gas dynamics case, so 





Vl v 2 

/iiA 1/2 wi jl 2 X l/2 v 2 







1 

-AV2d-i/2 y 





at) 


(41 

v . 





AV2d-v2 y 



a 1 


"2 









3± 



(3.42) 



det (/?+) and det ( " 2 

V P2 Pi 



are real nonzero quantities while oi^ are real, nonvanishing scalar functions 
of x. Note that block triangular structure plays a key role here as it allows 
the matrices to be block triangular. The RHS is equal to the product of 



(3.43) 
where 

(3.44) 



(3.45) 

and 
(3.46) 



det ( Vi I V 2 ) 












Vl 


v 2 







-fl 2 X 1/2 v 2 





* 


* 


1 


* 


* 


-X^ 2 d 



Vo 



\ 



( u 
-b^biu vi 





v 2 



\ 




* jliX l / 2 vi [i 2 X x l 2 v 2 

\ * * * x 1 ' 2 ^ 1 ' 2 ) 



det 



a+ 

or 
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Swapping columns, it's clear that the determinant (|3.43j) simplifies to 

/ u 

(3.47) — det vi i>2 —b^biu v\ V2 

det (-dV /2 (h) 1 ' 2 ) 

which is broken into kinematic and reaction pieces. The kinematic part then 
reduces as in Lemma LL2I to 



det (S+ ) det (vr W+ , eS+ ) det (vrW - ) , 

while the sign of the reaction part is (where we have included the minus sign 
from (|3~47jl ) 

— z 3 z 7 

whence the result follows. Nonvanishing follows from corresponding result 
for nonreacting gas dynamics and the fact the are nonvanishing for all 
x. " □ 

Nonvanishing of D(+oo) yields that the stability index T satisfies either 

f = sgn7 d A, 

or 

f = -sgn7 d A, 

as model parameters are modeled smoothly. Thus we find that the relative 
stability index, defined to be 

sgn7 d A, 

gives a measure of spectral flow. That is, changes in the sign of the relative 
stability index indicate a change in stability. 

In this case, however, we can do more and actually evaluate the absolute 
stability index, 

sgn J D'(0) J D(+oo), 

by relating our formula above for large A to the normalizations we've chosen 
at A = 0. 

Proposition 3.3. For A real and sufficiently large, 

(3.48) sgnD(A) = - sgn det(S + ) det(vrW + , eS + ) det(vrW")4^f U=o + 



Proof. The form of ()3.39j) allows us to connect to A = separately in kine- 
matic and reaction terms. Since the pair (A, B) satisfy the semidissipativity 
conditions, the result of the Appendix holds. As for the reaction term, pull 
back to A = is evident since "reaction" vectors have the form 
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hence projection onto the z component never vanishes, and thus cannot 
change sign. □ 

3.3. The Stability Index. Combining Propositions 13. 1 1 and 13.31 we obtain 
the result. 

Theorem 3.4. The stability index for a strong detonation with Lax 3-shock 
structure is 

f := - sgn-f d Aj NS Z3 Zj det(rf ,r^ : U^), 
where ^ms involves only the gas dynamics components. 

Proof. The stability index has the form T = sgn D'(0)D(+oo). From Propo- 
sition l3~T1 we find that D'(0) = 7^A. On the other hand from Proposition ^. 31 
the term det(7rW~) simplifies to det(rj~, r^ , U§) as in Appendix C, while 
the term det(7rW + , e§ + ) breaks into the product of 



7NS = det(^ 3 2 ) 

V*3,l S3.2/ 



and det(§ + ) the latter of which cancels with the factor det(S + ) appearing 
in Q3S). □ 

3.4. Reduction of A and Nonvanishing of A for Ideal Gas. The term 
A in the stability index has the form of a modified Lopatinski determinant 

A = det(rf ,rj, [U] + q). 

In place of the determinant above, it is more convenient for calculations to 
rewrite the quantity as the dot product of an appropriate left eigenvector 
and the jumps in the conserved quantities. Thus 

A = ls-([U] + Q) = l3-[U} + l3-Q- 
The left eigenvector l^ is 

( Pe , PeU 2 p e Pe 

In = \ p p e — cm H , c u, — 

V P 2 P P P 

while 

[u] = ([p],o,[£}) tr 

and q = (0,0, q) . It is then straightforward to calculate 

(3.49) A = (pp - — ) W + V -^f- (H + [!//>]) + —Me] + — ? 

V P J 2 P \P J p p 

Now using 

[ 4 +[i/p]=-— wii/p] 

p- p- 

in our expression for A and simplifying, we obtain 
(3-50) A = (-^[1/p] +P P -—) [P] + P - P+ [e] + P ~o 

\ v p j p p 

Similarly as in [HSj, we have 
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Claim 3.1. [e] + (p)[l/p] = — ^ where (p) = |(p+ + p_). 
Proo/. From (lRH3ll 

£] |"p 
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'£' 


+ m 


~p~ 


m 








.P. 




.P. 



or 



(3.51) 



[E] + 



u+q 
m 



Using the relationship E = \ + e, we can use (|3.51j) to calculate an expres- 
sion for [e], the jump in the specific internal energy. Thus we find 



(3.52) 



[e] 



m r -2l 



u+q 
m 



The term in parentheses in ()3.52|) can be simplified to 



(3.53) 



m r 



(P)[l/P]- 



□ 



The claim then follows by substitution into (|3.52l) . 

From Claim I3~T1 we find that 

— P+[e] = -—p+{p)[l/p] - —q, 
P P P 

and upon substituting into the expression for A, we find that the q terms 
cancel out. Therefore 



(3.54) 



[p] 



p e m' 



■[1/P]+PP- — + 3<P>) 
P P l J 



Since [p] 7^ 0, the condition A = may be written as 

cm p e f\p] 



(3.55) 
Next we note that 



Pp 



+ 

P P* 



M + (p) 



+ <p> 



p+ 







and thus Q3.55JI simplifies to 
(3.56) 



cm p e p + 

P P + — 3- 

P P 



0. 



Finally using the fact that the sound speed c satisfies c 2 = p p + p 2 pp e , we 
find that 

PeP+ 

= <r -t- 



. ■ . 2 1 Pe r i 

Pp H — = C + — [pj, 
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and thus we reduce (|3.56j) 

2 , Per l mc n 

c + — fa = 

p p 

P 2 P 

i P em2 n / l 771 n 
1 2~2" 1/P = 

p A c z pc 

Since the Mach number M satisfies M = — we obtain Majda's condition 
for inviscid shock instability 

(3.57) M 2 [l/p]p e - M -1 = 0. 

See |43| . We remark here that the jumps in the detonation formula above 
refer to endstates of the whole wave and not of the Neumann shock at the 
leading edge. Using that A reduces to the Majda condition for the stability 
of inviscid shocks, we now show that A does not vanish when the equation 
of state is assumed to be of the form 

P(P, e) = Tpe. 
We denote the compression ratio by 

P+ 

r = — , 
P- 

and we denote by r* the compression ratio of the Neumann shock. It follows 
immediately from the Rankine-Hugoniot diagram, Figure that 

(3.58) 1< r < r*. 

Also in the case of an ideal gas, we find similarly as 

(3.59) 1< r* < l±i = 1 + | 

7 — 1 1 

Combining Q3.58JI and (|3.59|) we find that 

2 

(3.60) < r - 1 < - 

Specializing the Majda condition to the ideal gas case yields, 

T(r - 1)M 2 - M - 1 < 2M 2 - M - 1 

= (2Af + 1)(M-1). 

The quantity (2M + 1)(M — 1) is nonpositive for — | < M < 1 and for 
strong detonations the Mach number M satisfies < M < 1. Finally since 
A = [p](r(r — 1)M 2 — M — 1) by our reduction above and since [p] < 
for 3-shock detonations, we conclude that A > for all strong detonation 
waves under the ideal gas assumption. 
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3.5. Evaluation of f in the Ideal Gas ZND limit. Here we use the 
structure of the singular manifolds contructed in the existence argument of 
|18j to determine the sign of the transversality coefficient 7^ and the other 
terms in the stability index. Recall that this coefficient is defined by the 
determinant 



Id 



det 



($1 

VCti 



(2,2 

^3,2 

C4,2 

c+' 

M,2 



C2,3 
^3,3 

(I* 



^2,6 
^3,6 
^4,6 
^4,6 



\ 



/ 



In order to use the structure of the singular manifolds, we must translate 
from the (^-coordinates of our stability analysis to the uTY Z-coordinates of 
the existence argument. This is accomplished in two steps. First, using the 
original C-coordinate transformation (|3.11|) and the A = eigenvalue equa- 
tion, we can connect the ^-coordinates to a set of intermediate coordinates: 
Pj, 8j, Zj, z'j via the linear transformation (dropping ±) 



(3.61) 



/C2A 


/ 621 












fpj\ 




hi 


633 












W 


d 





d 









\d' + &21 la 23d 


&2] 1 a 23d 





d) 







Note that the determinant of the linear transformation above is simply 



vm 



•c" 1 ^ > 0. 



hi hz d 

To obtain (J3.61|) note that from (|3.11f) we know 

(3.62) (2,3 = hiPj 

(3.63) = hipj + b 3 sSj 

(3.64) = dpj + dzj. 
Also from 1)3.64(1 it follows that 

Caj = d' pj + dp'j + dz'j . 

We take advantage of the second eigenvalue equation at A = 0, 

(a 2 ipj + a 2 ^£j)' = {hip'j)', 

and the fact that we are interested in the behavior of fast modes, so we can 
integrate up once to obtain 

Pj = b^lanPj + fc^a^j, 

which we can use to write p'j in terms of pj and Sj. Substituting into the 
equation for (' A ■ , we get 

(3.65) Ci,j = {d' + b2iOt2id)pj + b^l^dEj + dz'j. 
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Combining equations (|3.62j) - (|3.64j) and (j3.65|) yields (|3.61|) . 

The second step is to connect pj , £j , Zj , z'j , the variations in the conserved 
quantities, to Uj,Tj, Yj, Zj, the coordinates in which the construction of the 
singular manifolds has been accomplished. We use the relationships 

(3.66) pu = m 

(3.67) p(u 2 /2 + c v T) = 8 

(3.68) pY = z, 

and we linearize about the profile to obtain (note that subscript j indicates 
a variation while a bar indicates that a quantity is evaluated on the profile) 

(3.69) " - 
(3.70) 



Pj 
Si 



mu 



pu Uj 

c v pTj — pu~ 1 Uj(u 2 /2 + c v T) 



(3.71) zj = -pu^Yuj + pYj. 

Also since z' = p'Y + pY' , we linearize and use known relationships among 
the variations to obtain 

p'Yj + Y{b^la2\Pj + b^a^Ej) + Y'pj + p{- 



(3.72) z'j 
where Z is defined as in |18j by 

Z 



m m 
Td 3 W 3) ' 



Y- P ^ 
m 



Finally we can use (|3.69|) - ()3.71|) and (|3.72|) to write down the coordinate 
change as 



(3.73) 



(Pi\ 


/ 


—pu 1 








^ 




(Uj\ 




m 


- pu- l (u 2 /2 + c v T) 


c v p 
















-pu^Y 





P 









w 


\ 


Mi 


M 2 


P' + l 


m 
d) 







where 
(3.74) 
(3.75) 



M x = -pu-\b^a 2l Y + Y') + [m 
M 2 = b^laizYcvp 



pu 



u 2 /2 + c v T))b 2 ~ 1 1 a 2 3Y 



and we have dropped the bars on all terms in the matrix. Moreover, we note 
that the determinant of the matrix of this second coordinate change is 

c v p 6 u i — - > 0. 
a 

Therefore we find, by virtue of the fact that the determinants of the two 
change of coordinates matrices are positive, that 



(3.76) 



sgn 7d = sgn det 



fuf ■ 


■ u 6 \ 


r+ • 


■ T 6 


y? ■ 


■ Y 6 ~ 


W ■ 
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By a simpler calculation proceeding as above, we also find that 

(3.77) sgn 7ns = sgn det (p = sgn det (f + f + 

Finally the sign of T can be determined by a careful examination of the 
structure of the singular manifolds from which the solution is constructed. 

Recall from Chapter 1 that the structure of the singular manifolds is as 
in Figure 11.41 In this figure, the T-axis is perpendicular to the uZ-plane 
on the page, while the fourth missing direction in the phase space is the 
Y direction. Also the parabolic curve, C, on which the slow (reactive) flow 
takes place is not in a T = constant plane, see Figure 11.31 The diagonal 
dotted line indicates the intersection of the T = Ti plane and C. Thus the 
bold segments of the branches of C are below ignition temperature and there 
are no slow dynamics on those portions of the curve. Furthermore, we note 
that the unburned state (u + , T + , Y+, Z + ) is a degenerate rest point due to 
the ignition temperature assumption with a 3-dimensional stable manifold. 
On the other hand the burned state, (ii_,T_, 0,0), has a 2-dimensional 
unstable manifold (featuring a reactive and a kinematic direction) and a 
2-dimensional stable manifold. Note that the T and Y directions are both 
stable directions, see the discussion of existence in Chapter 1 particularly 
equations (fTT7|) - (fl~19|) and equations (|L25|) - (|L28|) . 

We evaluate 

(3.78) - sgn7d7 NS det(rf , r%, U e )z^z x \^ OQ 

at the "corner" where the fast manifold which approaches the burned state 
intersects with the opposite branch of C. In figure l3~Tl we see a schematic 
indicating the relevant vectors in the calculation. The 1, 7 arrow corresponds 
to the profile, while the 6 arrow corresponds to the kinematic direction at 
— oo. The dashed 2 arrow represents the stable manifold in the T-direction, 
while the curvy 3 line corresponds to the missing stable Y-direction. Then 




Figure 3.1. The Point of Evaluation 
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we find that sgn7d7NS can be computed as 

/+ -\ 

(3.79) sgndet 



<1<1 ( + ° ] I. 



- - - 
* - 

\+ o o o/ 

These choices force 

(3.80) sgndet(rf,rJ J l7 6 -) = -l J 

To see (|3.8Uj) . note that det(rj~, r^, Uq |_oo) = 1$ • Uq where 1% is the ap- 
propriate left-eigenvect or, (simplified due to ideal gas assumption) 

/ Tu 2 N tr 

(3.81) 1% = [-CU + — ,c-Tu,T 



Then since 

(3.82) U 6 = (p 6 ,m 6 ,S G ) tr = (p 6 ,0,£s) tr , 
we obtain by combining ([3. 81ft and (|3.82j) that 

Tu 2 

(3.83) / 3 - • = (-cu + — )pe + T£§ . 

To take advantage of the signs we know, we translate into u, T coordinates. 
We use (|3.61|) and (|3.73j) in (|3.83j) . and then some elementary simplification 
yields 

(3.84) Z 3 ~ ■ U Q = (cp - ^) n 6 " + Tc vP T 6 . 

To evaluate the sign of 1)3. 84|) . we note that because we consider a 3-shock, 
u < 0, and thus the coefficient 

/ CvP\ 

is positive. Finally since both sgn-Ug and sgnT 6 ~ are —1, (|3.8Uj) follows. 
Also we find then that 

(3.85) sgnz^ = sgnpY 3 + = —1. 
Lastly, we find that 

(3.86) sgnzf = sgn^l-oo = sgn(py a; |_ 00 - pvr l Yu x \-oo) = +1. 

Tracking the signs computed in (|3,79|) - (|3.86[) and combining with (|3.78|) and 

the fact that we computed 

A > 0, 

for all ideal gas strong detonations, we discover that in the ideal gas ZND 
limit, the stability index satisfies 

f = -sgn7rf7 NS Adet(r 1 ",r^,[/ 6 ~)z3 f z :E |_ 0O 
= (-1)(-1)(+1)(-1)(-1)(+1)=+1 

which is consistent with stability. This completes Theorem 11.51 
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Finally, we also remark here that under the ideal gas assumption we found 
sgnA by comparison to the Neumann Shock. Restating this, we found 
that if the Neumann Shock is "index stable," i.e. satisfies the stability 
index necessary criterion, then the corresponding strong detonation is also 
"index stable." A natural question then is: under what conditions on the 
equation of state does this remain true? Or perhaps more importantly, 
what is the actual stability relationship, not just the relationship between 
stability indexes, between the shock and the strong detonation? The partial 
information gathered from the stability index approach definitely motivates 
further investigation into this question. 

3.6. Stability Index for Multiple Reactants. In actual combustion, the 
chemical reactions involved are typically more complicated than a single one- 
step reaction. It's natural to model reactions in an s + 1 component gas by 
a system of the general form 

(3.87) U t + f(U) x = (B{U)U x ) x + Q<5>(U)z 

(3.88) z t + (v(U)z) x = {D l (U)z x ) x + {D 2 (U,z)U x ) x -^(U)z, 

where U = (p, m, £) tr and / and B are as in the Navier-Stokes model for 
gas dynamics considered above. The vector z = pY £ R s measures the 
quantities of each of the reactants; the constant matrix Q £ R 3xs records 
the heat released in each reaction, thus 

/ ••• 
Q= ••• 

\ <?i ••• q s 

The positive definite matrix 3>(J7) € R sxs incorporates the reaction rates; 
the diffusion matrices satisfy D 1 € R sxs and D 2 € R sx3 ; and the scalar 
function v(U) is simply ^ = u, the velocity. 

It turns out that the number of equations required to model the chemistry 
can be reduced through the use of progress variables. See ^1] or [3U] for 
further discussion. Our analysis largely applies to the more general multiple 
reactant case with the following important exception below. 

In the more complicated multi-species case the large A calculation for 
the reaction block is not so straightforward. In particular the connection 
to A = is not clear; we indicate here some partial results along this line. 
Examining the reaction portion of the characteristic equation on the +oo 
side, we find since the reaction function satisfies $ = 0, 

(3.89) (f?D - pV + \I)z = 0. 

Since the convection in the reaction equation is just the background velocity, 
V is scalar, hence it commutes with the species diffusion. This implies that 
there is a symmetrizer V° such that V°V and V°D are symmetric, and 
Re V°D > 0. Thus arguments using Theorem IA.4I apply. On the other 
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hand, things are trickier on the — oo side the characteristic equation reads 

(3.90) (fj?D - pV + XI - $)z = 0. 

Here, we point out that the argument of [3], Lemma 7.2 goes through word 
for word for constant-coefficient eqautions of the form Xz + Az' + Bz" + Cz = 
provided that A is symmetric and B and C are positive definite, i.e. 
Re-B,ReC > 0, (more generally there is a coordinate transformation for 
which this is true) with the slightly modified Lyapumov function H(v, uu) = 
(l/2)(Aw,w) + Re(x, (A + C)v). Indeed the conclusions become somewhat 
stronger extending to all Re A > and not only Re A > as in In our 
case A = v is scalar, B = D and C = $, so we find that the machinery of 
[2] may be applied provided that there exist coordinates in which D and $ 
are both positive: in particular if D is scalar. 

3.7. Small q Stability. Here we examine the q — * limit, and take advan- 
tage of the simplification in the equations when q = 0. Using a continuity 
argument, we prove strong spectral stability, that is, nonexistence of eigen- 
values with Re A > and A ^ 0, as well as transversality of the connecting 
profile jd / 0, and low-frequency stability A / 0. It is expected that by 
following the program of 68 and jlBlETJEH]) these three properties should 
be sufficient to conclude full nonlinear orbital stability. Such a result would 
be an extension to the reacting Navier-Stokes model of the results of 3H].j32j 
in which the authors show nonlinear stability of strong detonations in the 
Majda model as q —* 0. We note that in both [3E| and |32j . as here, working 
with the integrated equations is a key ingredient in the analysis. 

We rewrite the system (|l.l|) - (ll.4|) using U to represent the vector of gas 
dynamical variables and z = pY in the first three equations, 

(3.91) (U + qz) t + (/([/) + quz) x = (B(U)U X ) X , 

(3.92) (pY) t + (puY) x = (pdY x ) x - k<p(T)pY. 

Now we define U by U = f^^iU+qz), and consider the linearized eigenvalue 
equation from Q3.91JI in the integrated variable U, 

(3.93) XU + AU' -{A- u)qz = BU" , 
and from (|3~92|) 

(3.94) X{pY + pY) + (puY + puY + puY)' = {palY 1 + pdY x )' 

— k(ppY — k(ppY — k(p (T)TpY. 

Note that when q = 0, ()3.93j) reduces to the integrated eigenvalue equation 
for gas dynamics about a gas dynamical profile. 

Lemma 3.5. The limiting system (in integrated form) has no eigenvalues 
with Re A > provided the limiting shock is stable. (In particular, the lim- 
iting shock is stable if the amplitude is sufficiently small.) 
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Proof. Since the limiting shock is spectrally stable, we have by standard 
considerations [HE] that the integrated eigenvalue equation for gas dynamics 
supports no eigenvalues on Re A > 0, and so (p, u, £) identically vanish for 
any eigenfunction of the limiting eigenvalue equations as q — ► 0. Thus 1)3.94(1 
has the simpler form 

X(pY) + (puY)' = (polY')' - kCppY. 

We note that pu = m is real and constant on the profile. Then taking the 
standard complex I? innerproduct of the above equation with Y, we have 

(X(pY),Y) + ((mY)',Y) = ((pdY')',Y) - (k<ppY,Y), 

so that integrating by parts, and taking real parts yields 

(3.95) Re X(pY, Y) + + Re{k<ppY, Y) + Re{Y f , pdY') = 0. 

since (mY', Y) = {(mY)' , Y) = -(mY, Y') = (mY, Y') so that ((mY)', Y) is 
purely imaginary. But for 1)3.95)) to hold, we must have either Re A < or 
Re A = and also Y' = which implies that Y is constant. But then we 
need also Ke(k(ppY, Y) = 0, so that the constant value for Y must be 0. We 
conclude that the only nontrivial solutions must correspond to Re A < 0. 

In the case that species diffusion is neglected things are even simpler; we 
write the reaction equation as 

(3.96) \z + (vz)' = -cpz, 

where z = pY as usual. Applying the Gap Lemma to (|3.96f) . we find that 
behavior at +oo is governed by the limiting constant-coefficient equations. 
These are easily seen to support no stable modes. Indeed, at +oo, we have 
(p = 0, so there is no reaction at all. Then the equation 

\z = —(vz)', 

which can be rewritten in terms of the new variable w = vz as 
(3.97) 

is obviously blowing up for Re A > 0. □ 

Corollary 3.6. For sufficiently small q, the integrated eigenvalue equations 
associated with ()3.91|) and (|3.92(1 have no Re A > eigenvalues provided that 
the limiting shock is stable. 

Proof. The respective Evans functions vary continuously, and the limiting 
Evans function is nonvanishing on Re A > 0. □ 

Proposition 3.7. // the limiting shock is stable, then then small-q detona- 
tions are strongly spectrally stable. That is, they have no eigenvalues for 
Re A > and A ^ 0. Moreover, there holds D'(0) = 7^A ^ 0, where 7^ and 
A are as defined in (|3.26f) and (|3.25|) 
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Proof. When A / 0, we may integrate the divergence form gas eigenvalue 
equation to deduce that (U + qz) has zero integral, and thus U, defined 
as (U + qqz) lies in L? if U,z do since by the gap lemma functions 
decay exponentially if at all as do their integrals [HE]- Thus existence of 
an eigenfunction for the linearized eigenvalue equation and its integrated 
version are equivalent. Existence of a transverse connection 7d 7^ follows 
likewise by continuity from the result for the limiting equations, provided 
there exists a transverse connection for the limiting gas-dynamical shock. 
Finally A/0 follows by inspection from the corresponding property for the 
limiting gas-dynamical shock, since A for q = reduces to this case, which 
property is a necessary condition for stability [SHI- ^ 

The results of this section clearly extend to the case of multiple reactants. 
Indeed for the zero species diffusion case, the number and type of reactants 
plays no role. For the D 7^ case, the arguments above carry through if D 
and $ are simultaneously positive, or if there exists a constant coordinate 
change making them both positive. In particular the argument applies if 
the diffusion D is scalar and $ = ipK for K constant and ip scalar. 
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Appendix A. Real Viscosity 

For completeness and convenience, we provide here a revised version of an 
appendix of [US] using 6(>] and including a discussion of the extreme shock 
case in Section fA,2l below. Systems modeling gas dynamics have the general 
form 

(A.l) U t + F{U) X = {B(U)U X ) X , 

where 



v J ' \ 9 r \hb 



>2 



and 



uJeR n ~ r , v,g,eR r , iieR rx (">, 6 2 eR m . 

We note that in the isentropic gas dynamics case n = 2 and r = 1 while in 
the Navier-Stokes case n = 3 and r = 2. Our interest is in traveling wave 
solutions of the form 

(A.2) U = U{x), lim U(x) = U± = (u±, v±). 

x— >±oo 

Standard assumptions for equations in this generality are: 

HO: F,B £ C 2 . 

{(i)Refj(6 2 ) > and 
(ii)(^) full rank on {U(-)}, moreover 
(iii)(^fu - foi(6 2 ) _1 d/ 1) ) is real. 
H2: a (dF(U±)£) real for £ G R, a(dF(U±)). 
H3: Re a (£dF(tf±) - ^B(U ± )) < for £ G R. 

H4: Solutions of (|A.l[) - l|A.2j) form a smooth manifold {u 5 }, 5 G W C R f . 

These hypotheses are analogous to those of the strictly parabolic case con- 
sidered in e.g 65 a , with (Hl)(i) and (HI) (iii) ensuring local well-posedness. 
Indeed, they are the standard set of conditions identified by Kawashima 30 ; 
for further discussion, see [HS1- The condition (Hl)(ii), may be motivated 
by consideration of the traveling wave ODE 

(A.3) /(«,«) = /(«_,«_), 



(A. 4) b\u + biv = g(u, v) — g(u-, 

For, (Hl)(ii) is readily seen to be the condition that (jA.4j) describes a non- 
degenerate ODE on the r-dimensional manifold described by (|A.3|) : thus, 
this is a reasonable nondegeneracy condition to impose in the study of vis- 
cous profiles. Condition (HI) (iii) also arises in the analysis of the eigenvalue 
ODE; see the discussion of consistent splitting in Appendix A2, [HE]- In the 
symmetrizable case it holds automatically. 

We remark, finally, that (Hl)(ii) (indeed, all of hypothesis (HI)) is satis- 
fied for gas and plasma dynamics precisely when particle and shock velocities 
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are distinct, which is always the case along a shock; for a study of viscous 
profiles in these contexts, see [2^, (20], [T6| . 

Let i + denote the dimension of the stable subspace of i_ denote 

the dimension of the unstable subspace of df l {u^), and i := i + + i_. Let 
eL|_ denote the dimension within the submanifold / = constant of the stable 
manifold at (u+,v+) of traveling wave ODE (JA.4JI . and d- the dimension of 
the unstable manifold at (u_,i>_), and d := d- + Then, we have the 
following result analogous to that of Majda and Pego [MP] in the strictly 
parabolic case. 

Lemma A.l. Under assumptions (H0)-(H3), (u±,v±) are hyperbolic rest 
points of the reduced traveling wave ODE ()A.4|) . In particular, traveling 
wave solutions satisfy 

(A.5) \(d/dx) k ((u(x),v(x)) - (w±,u±))| < Ce~ 9 ^, fc = 0,...,4, 

as x — > ±oo. Moreover, the type of the connection agrees with the (hyper- 
bolic) type of the shock, in the sense that 

(A. 6) d — r = i — n. 

Proof. Integrating (|A.4|) from — oo to x and rearranging, we may write 
(|A.3|) - (|A,4|) in the alternative form 

Linearizing (|A.7|) about U±, we obtain 



(A.8) 

or, setting 



the pair of equations 



and 



u\ fu fv\ f \ /u 



z i\ ._ ( fu fv\ ( u 



Z2j ' \bi b 2 J l{u±) \v 



(A.9) A = {9u 9v )[il £)" {i) m ^> 

the latter of which evidently describes the linearized ODE on manifold (jA 
Observing that 
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by (H2) and (Hl)(i), we find that the coefficient matrix of (|A.9|) has no zero 
eigenvalues. On the other hand, it can have no nonzero purely imaginarly 
eigenvalues i£, since otherwise 



fu fv 

9u 9v 



and thus 



fu fv 

9u 9v 



fu 

h 





h 



fv 

h 



o 



fu fv 

b\ b 2 



-i 



for ^ / G R, in violation of (H3). Thus, we find that U± are hyperbolic rest 
points, from which (|A.5|) follows. Relation (jA.6|) now follows from Lemma 
lA~5l below. □ 



The linearized eigenvalue equations about U(-) are: 
(A.10) 
and 



(A.ll) 
where 



(Anu + Auv)' = -An, 
in' + b 2 v'Y = (A 2 \u + A 22 v)\ 




h 





h 



B(U), 



An 
A 2 i 



A 12 
A 22 



U := dF{U)U - dB(U)(U, U' x ] 



and 



denotes d/dx; in particular, note that 
(A 11 ,A 12 ) = df(U). 
Utilizing the change of variables 



(A.12) 



z 2 



An 



A 12 
b 2 



we can write the eigenvalue equation as a first order system 

(A.13) Z' = k(x,X)Z, Z = ( Zl ,z 2 ,z' 2 ) tr . 

We note here that z 2 has r components. 

The consistent splitting hypothesis can be verified by a limiting analysis 
as A — > +oo, carried out without loss of generality in original coordinates 
W, for which the asymptotic characteristic equations become: 

Mil + A 

ll2 



(A. 14) det 
This yields n — r roots n 
(A.15) 



Ml2 
.21 



jiAi 2 - jj?b\ nA 22 - [i' 2 b 2 + A y d 



jlX, [~l ~ 1, satisfying 



Mn + 1 

h 



A 12 
b 2 
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or 

(A.16) -/T 1 G a (A n - Anb^b^, 

and 2r roots [i ~ /2A 1 / 2 , /i ~ 1, satisfying 

^ U -aw) ± (:H). 

or 

(A.18) /T 2 G a(6 2 ). 

By assumption (Hl)(iii), ()A.16|) yields a fixed number k/(n — r — k) of 
stable/unstable roots, independent of x, and thus of ±. Likewise, (Hl)(i) 
implies that (jA.18|) yields r stable, r unstable roots. Combining, we find the 
desired consistent splitting, with (k + r)/ (n — /c) stable/unstable roots at 
both ±00. We can thus define an Evans function as usual as 

(A.19) D(X) = det(Z+ . . . Z+ +r , Z~ +r+1 , . . . Z^ r \ x=Q 

where {Z^ , . . . Z^,}, {Zj~ +r+1 , . . . Z~ +r } span the stable manifold at +00, 
unstable manifold at —00 of (|A.13|) . Notice that the Evans function in Z 
coordinates is just a constant multiple of the corresponding Evans function 
defined in W = {u,v,v') 1 coordinates. 

A.l. Stability Index. The stability index is defined to be 
(A.20) f := sgn(d A ) £ £>(0)£>(+oc). 

The low frequency calculations of D'(0) are detailed for the Navier-Stokes 
model in Section 3. Note that I = 1 in this case. Here we evaluate the sign 
of D (A) as A — > +00 along the real axis. 

Lemma A. 2. Let D{-) denote the alternative Evans function 
(A.21) D(X) := det(W+, . . . W+ +r , Wj^. +V . . . W~ +r )\ x =o 

computed in the original coordinates W . Then, for real A sufficiently large, 
there holds 

(A.22) sgnZ)(A) = sgndet(S + , U + ) det(vrW + , eS + ) det(eU~, ttW~) + 0, 

where ir denotes projection of W = (u,v,v') onto (u,v) components, and 
S(x), U(x) are real bases of the stable/unstable subspaces of (An — A^b^bi) 
(note: (n — r) dimensional), with eu := (u, — b^biu) denoting extension. 

Proof. Recalling (Hl)(ii), we have that the coordinate change(u, v) — > (z\ , z%) 
is invertible, and so we may work equivalently in (u, v, 23) coordinates, where 
23 := b\u' ' + b2v' '. Then, we find from l|A.15|) . ()A.17|) that the stable/unstable 
manifolds of the frozen eigenvalue equation at any (fixed) x are spanned by 
vectors of form 

u \ 
-b 2 hu , 
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with u an unstable/stable eigenvector of (An — A^b^bi), — the corre- 
sponding eigenvalue; and vectors 



v 

with v an eigenvector of b2, —fj-~ 2 the corresponding eigenvalue. Rescaling 
and applying the tracking lemma, we thus obtain 



(A.23) D(X) ~ det 



u 

-b^biu 



■•• 
v ■ ■ ■ v 

-fi\ l ' 2 v ■■■ fi\ l ' 2 v 



2r 



x det 








2 

a + 








al 





2 

at 


& 





Pi 








0L 





Pi 



where both 



det 



on, a 2 + 

Pi Pi 



and det 



al a 2 _ 
P 2 Pi 



are real, nonzero quantities. The right hand side of (|A.23|) can be rewritten 
as 



(A.24) 



det(S + ,U + )det(Vrdet J 



1 2 

a, al 



Pi Pi 



det 



al a 2 



/3 2 0L 



i » 



where by det(V) we refer to the r x r determinant coming from the u com- 
ponent of (|A.23|) . On the other hand, the term 



det(vrW + ,eS H 



in HA.22JI can be simplified to 
(A.25) det(7rW + ,vrU^ 





a\ 





Pi 


Pi 














where (W+,U+) = W+, but det(vrW + , vrU + , eS + ) has the form 

U+ S+ 



det 





V 



Finally we see that (jA.25|) is simply 
(A.26) det(V) det(S+, U+) det 

Similarly, the term 

detain, vrVT 



ai a 2 + 
Pi Pi)' 
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simplifies to 

fa 1 a 2 

(A.27) det(V)det(S~,U~)det f ^ qT 

combining l)A.26|) and (|A.27|) we find since sgn det(S + , U + ) = sgn det(§~, U~) 
that the expressions (|A.22f) and (|A.23|) agree modulo the real, positive factor 
det(S+,U+) 2 . □ 

We further make the assumptions of semidissipativity 

(+) There exist symmetrizers A± such that A-t- are symmetric and 
ReA° ± B± < 0, 

and block structure 



(++) (Aiyi'B ± {A 



± 



'0 

Both of these assumptions hold for the compressible Navier-Stokes equa- 
tions. When they hold, more can be said about the sign of -D(A) for large 
real A . 

Lemma A. 3. Assuming (+)-(++), there holds, for sufficiently large, real 
X: 

(A.28) 

sgnJD(A) = sgndet(S + ,U + )det(7rW + ,eS + )det(eU",7rW") u=0 ^ 0, 

where D(-) as in \A.21\) denotes the Evans function computed in original 
coordinates W = (u,v,v'y. 

Proof. Without loss of generality, we may take A symmetric, B = ( _ , 

and Re 62 > 0, by the transformation 

A^(A°)7A(A°)-v, B -» (A°)*B(A )-?. 

It is sufficient to show that quantity ()A.28|) does not vanish in the class (+)- 
(++). For, since D(X) does not vanish either, for real A sufficiently large, we 
can then establish the result by homotopy of the symmetric matrix A± to an 
invertible diagonal matrix (straightforward, using the unitary decomposition 
A = UDU*, U*U = I, and the fact that unitary matrices are homotopic 

either to I or —I) and of B± to ^ j ^j (e.g., by linear interpolation of 

the positive definite b 11 to I r ), in which case it can be seen by explicit 
computation that (|A.28|) is independent of A G [0, +00]. We note that the 
endpoint of this homotopy is on the boundary of but not in the Kawashima 
class, since eigenvectors of A are in kernel of B; indeed, our definition of 
semidissipativity is is not the "strict" dissipativity condition of Kawashima, 
but a nonstrict version. However, it suffices for the present, purely linear- 
algebraic purpose. 

We begin by examining det(-7rW + , eS + ). When A = 0, a bifurcation anal- 
ysis as in of the limiting constant-coefficient equations at ±00 shows 
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that the projections ir of slow modes of W + may be chosen as the unstable 
eigenvectors of A, corresponding to outgoing characteristic modes, and 
the projections of fast modes as the stable (i.e. Re \i < 0) solutions of 



vj \0 

or without loss of generality 

'A u A12 \ fu\ fO 



and thus of form 



An A 2 2 - fJ.b 2 J ± \v J \0 
(An^Auv 



v 



where 



{b 2 1 {A 22 - A 21 {A u y 1 A 12 ) -nl)v = 0. 

Likewise, using b\ = 0, we find from the definitions of §, e in the statement 
of Lemma IA.2I that stable solutions S + are in the stable subspace of An , 
with e S + = ( S g ) , hence vectors eS + lie in the intersection of the stable 
subspace of A and the kernel of B. Our claim is that these three subspaces 
are independent, spanning C n . Rewording this assumption, we are claiming 
that the stable subspace of (A)^ 1 B, the center subspace ker B intersected 
with the stable subspace of A, and the unstable subspace of A are mutually 
independent. (Note: that dimensions are correct follows by consistent split- 
ting). But, this follows by Lemma lA.41 below. Similar considerations apply 
to det(eU _ ,7rW-). □ 

A key to the calculations above is the following modified lemma of Serre 
from [HE]- We begin by fixing notation. For a matrix M, we denote by S(M) 
and U(M) the stable and unstable subspaces of M. Also we we have 

M V {M) = {v g V | (v,M) < 0} 

and 

V(M) = {v | (v,Mv) > 0}. 

Lemma A. 4 (Modified Serre's Lemma (HSD- Let A be a symmetric, invert- 
ible matrix and let B be a positive semidefinite matrix, Ke(B) > 0. Then 

(1) the subspaces S{A~ l B) @ Nh erB (A) andU(A) are transverse. 

(2) the subspaces S{B A -1 ) @ Nh erB (A) andU(A) are transverse. 

Proof. It suffices to prove the first claim, as the second follows immediately 
from the first due to the similarity transform B — > ABA^ 1 . Since A is 
symmetric, the unstable subspace U(A) is equal to the subspace V{A), which 
for notational convenience we denote by V . Suppose, in order to obtain a 
contradiction, that xq ^ lies both in the subspace S{A~ 1 B) (B M^ erB (A) 
and in V. Thus we suppose that 

Xq = X\ + X 2 
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where x\ G S(A 1 B) and %2 G A/jj erB (A) . Define by the ordinary 
differential equation 

(A.29) x = A~ l Bx, x(0) = x . 

It follows then that x(t) — > xi as t — > +00 and thus 
(A. 30) lim (x{t),Ax(t)) < 0. 

However, 

(x,Ax)' = 2(A- 1 Bx,Ax) 
= 2(Bx,x) 
> 0. 

This implies that (xq,Axq) < which is a contradiction to our assumption 
that xq belongs to V. □ 

A. 2. Evaluation of T. Using the z-coordinates we may regard the traveling 
wave ODE as an r dimensional first order dynamical system. We denote by 
d± the dimensions of the stable manifold at Z2+ and the unstable manifold 
at Z2- ■ We also define d to be the sum d = d+ + c?_. It follows then that 
1 < d± < r. Also we denote by i± the number of characteristics entering the 
shock from the left ( — ) and the right (+). We put i = Corresponding 
to [H], we have: 

Lemma A. 5. The following relations hold 

(1) n — i + = r — d + + dim^/(^4 + ) 

(2) n - i_ = r - cL + dimS(A^) 

where A = (An — A^b^bi) . Moreover n — i = r — d. 
Proof. Equating the dimensions of Z + at A = and as A — > 00 we find 
d\mU(A + ) + d + = dimU(A + ) + r, 

or 

(n — i+) + d + = dimU(A + ) + r. 
Similarly as x — > —00, we find 

(n — z_) + d- = dim5(A_) + r. 

That n — i = r — d follows by adding the two equations and noting that 
dimU(A+) and dimS(A^) are constant and sum to n — r. □ 

Corollary A. 6. For (right) extreme shocks, i+ = n, there holds also d+ = r. 
Thus the connection is also extreme and d\m.lA(A) = 0. 

Proof. This follows at once from Lemma IA.5I due to the fact that d+ < r 
and dimU(A) > 0. □ 
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The import of Lemma IA.5I is that the "parabolic" and "hyperbolic" types 
of connections agree. From Corollary IA.61 we may deduce that 7 for an 
extreme right (i.e. n-shock) Lax profile consists of a Wronskian involving 
only modes from th +00 side, and is therefore explicitly evaluable. For, 
working now in (zi,Z2,z' 2 ) coordinates, we obtain 7 as a determinant of z% 
components only. 

Moreover, the expression (|A.28|) simplifies greatly. In the {z\,z%) coordi- 
nates we have that U = 0, S is full dimension n — r, and eS consists of vectors 
of the simple form (z,0). This means that det(S, U) simplifies to just detS, 
while det(-7rZ + , eS + ) simplifies to the product of detS + and 7. Therefore, 
this term, similarly as in the strictly parabolic case, cancels with term 7 in 
the computation of the stability index. 

Finally, det(eU~, 7rZ~) simplifies to det(rj~, . . . , r~_ 1? u'), times the de- 
terminant of the coordinate transformation from W to Z coordinates, the 
latter determinant cancelling with a like factor appearing above. We are left 
in the end with the following very simple formula. 

Proposition A. 7. In the case of an extreme right (n-shock) Lax profile, 

(A.31) f = sgndet(rr, . . .,r~^, [«]) det^, . . . , r"_ l5 u' /\u'\(-oo)) 

We emphasize that this is identical with the stability index in the strictly 
parabolic case. The only very weak information required from the connection 
problem is the orientation of v! as x — > —00, i.e. the direction in which the 
profile leaves along the one-dimensional unstable manifold. We remark that 
in the case of isentropic gas dynamics, the traveling wave ODE is scalar, and 
thus the orientation of v! is determined by the direction of the connection. 
See Section 2.4 or IIU] for further details. 
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